library(ggplot2)
library(pander)
library(knitr)
library(tidyr)
library(mgcv)
library(diversitree)
library(repmis)
library(lme4)
library(ape)
library(geiger)
library(grid)
library(gridExtra)
library(nlme)
library(phytools)
library(brms)
library(ggridges)
library(caper)
library(purrr)
library(reshape2)
library(ggExtra)
library(grid)
library(dplyr)
library(BAMMtools)
library(gridExtra)
library(car)
library(coda)
library(MuMIn)
library(parallel)
library(HDInterval)
library(phangorn)
library(kableExtra)
library(RColorBrewer)
library(tibble)
library(stargazer)
library(ggraph)
library(stringr)
#devtools::install_github("Ax3man/phylopath")
library(phylopath)
library("EBImage") # for images
library("ggtree")
source("bamm_extinction/functions/check_and_fix_ultrametric.R")
source("functions/essim.R")
source("functions/diversification_rate_calculation.R")
source("functions/Vdodge_function.R")
source('functions/symmetrise_scale.R')
source('functions/Correlation_matrix.R')
source('functions/gheatmap2.R') #Gets rid of spaces between tiles
knitr::opts_chunk$set(cache=TRUE)

Compiling the data

We aimed to compile data on Passerine birds, which represent the largest group for which relatively complete phylogenetic, trait and environmental data are available. The following details the source and use of the data.

Speciation Rate

To obtain speciation rate estimates we used Bayesian Analysis of Macroevolutionary Mixtures (BAMM). Additionally, we compared our BAMM output files (event data) to an existing rate dataset (Harvey et al. 2017).

Rather than obtaining one estimate of speciation rate from one tree, we ran BAMM 100 times on 100 trees plus an MCC tree to obtain uncertainty estimates of the speciation rate generated from the phylogeny produced by Jetz et al. (2012). We also estimated extinction rate from BAMM and speciation rates using the diversification rate and node density statistic (\(\lambda_{DR}\) and (\(\lambda_{ND}\)). See Title and Rabosky (2018) for a comparison of methods.

Sexual dichromatism and Male-biased Sexual Selection

The first proxy for sexual selection that we used was based on measures of sexual dichromatism. We used two existing datasets; a complete set of male and female plumage scores in Passerines from Dale et al. (2015) and reflectance data from 1000 birds from Armenta et al. (2008).

Briefly, Dale et al. (2015) obtained male and female plumage score, the ratio of which provides an index of sexual dichromatism (SDi). The plumage score was extracted using the Handbook of the Birds of the World Del Hoyo et al. (2011). By scanning images of males and females in this book across multiple patches, Dale et al. (2015) obtained mean plumage scores from RGB values. Dale et al. (2015) compared their RGB scores to a dataset of reflectance measurements of Australian birds. The two datasets were correlated but the relationship was expectedly noisy (R2 = 0.67). Here we compare the Dale et al. (2015) data against another reflectance dataset from Armenta et al. (2008), which estimated sexual dichromatism as the mean number of just noticeable differences (JNDs) between male and female plumage colour (discriminability), based on a model of bird colour vision.

Alongside sexual dichromatism, Dale et al. (2015) compiled available data for traits such as: Body size, tropical life history, Sexual selection, Cooperative breeding and Migration. As a second measure of sexual selection, we used a phylogenetic PCA (PC1) that estimates male-biased sexual selection from a species level dataset of sexual size dimorphism, social polygyny and paternal care. Higher values indicate higher size dimorphism (larger males), higher polygyny and less paternal care.

Environmental predictors

We obtained species range distribution maps from Birdlife International (BirdLife International and Handbook of the Birds of the World 2017). These species range maps cover nearly all species of birds. From these spatial information files we obtained estimates of the following (note that not all the data extracted was subsequently used in the analysis):

  • Range Size
  • Average and standard deviation in 19 bioclimatic variables (each range randomly sampled 1000 times)
  • Average and standard deviation in 19 bioclimatic variables from the last-glacial maximum (LGM) and the last-inter-glacial (LIG)
  • Net primary productivity (NPP) estimates and variability
  • Average and standard deviation of human population density in the species ranges.

Below we outline the code used for extracting spatial data from bird ranges, however no raw data is provided alongside this file due to the large file size of the shapefile and raster information. We have put the extracted environmental variables in a csv filed called complete.dataframe.csv for convenience. Briefly for each of the ~ 6,000 species, we:

1) Randomly sampled the range polygon 1000 times:

#Increse the iterations (defaukt is 4) so we can obtain complete samples of each range
bird.points <- lapply(bird.ranges, 
                       function(x) {spsample(x, n=1000, type="random", iter = 30)})

2) Extracted the bioclim or other (e.g. NPP) data from each of the points:

bird.values <- lapply(bird.points, function(x) {raster::extract(bioclim_data, x)})
bird.values <- as.list(data.frame((bird.values)))

3) Summarised the extracted data across the 1,000 points into a summary value of interest (means and sd) for each species and exported that data.

#Obtain means per variable per species
bird.frame <- lapply(bird.values, function(x) {as.data.frame(x)})
bird.summary <- lapply(bird.frame, function(x) {
(as.data.frame(apply(x,2,mean, na.rm =T)))})

#transpose
bird.means <- t(as.data.frame(bird.summary))
bird.means <- (split(bird.means, 1:19))

#Now add column of species name: Same order carried through
bird.means <- cbind.data.frame(shps.jetz, bird.means)

rownames(bird.means) <- NULL
colnames(bird.means) <- c("binomial","bioclim1", "bioclim2", "bioclim3", "bioclim4", "bioclim5", "bioclim6", "bioclim7", "bioclim8", "bioclim9", "bioclim10", "bioclim11", "bioclim12", "bioclim13", "bioclim14", "bioclim15", "bioclim16", "bioclim17", "bioclim18", "bioclim19")

#Write csv
write.csv(bird.means, 'data/bird.means.csv')

4) The process above was repeated when we used gridded data for the LGM, LIG, NPP and human population density, while range size was obtained from the following code:

bird.range.size <- sapply(bird.ranges,
                       function(x) {(area(x))})
bird.range.size <- sapply(bird.range.size, function(x) {sum(x)})
bird.range.size <- as.data.frame(bird.range.size)
bird.range.size <- cbind.data.frame(shps.jetz, bird.range.size)
rownames(bird.range.size) <- NULL
colnames(bird.range.size) <- c("binomial", "range.size.m2")

#write.csv
write.csv(bird.range.size, 'data/bird.range.size.csv')

Generating biologically relevant predictors

In this study we assess the relationship between sexual selection and extinction risk. In doing so we attempt to take into account as many other predictors of extinction as possible, primarily through environmental variables. Our model structure seeks to contain:

Extinction/Diversification ~ Sexual selection
+ Range size
+ Short temporal variability of temperature (mean BIOCLIM4)
+ Spatial variability of temperature (PCA1) [residual.PCA1]
+ Long-term variability of temperature (LIG)
+ NPP

These predictors can be obtained from the compiled datasets seen here:

plumage.scores <- read.csv('data/plumage_scores.csv')
#Generate sexual dichromatism score: 
plumage.scores$SDi <- abs(plumage.scores$Male_plumage_score - plumage.scores$Female_plumage_score)

#Make DF with enviro variables and plumage scores
complete.dataframe <- read.csv('data/complete.dataframe.csv')
complete.dataframe <- left_join(plumage.scores %>% dplyr::select(binomial, Male_plumage_score, Female_plumage_score, SDi), complete.dataframe, by = "binomial")

Environmental Spatial Variability PCAs

To obtain estimates of spatial variability in environment we performed PCAs of the standard errors of the bioclimatic variables across the 1000 random points extracted per species (excluding seasonality in temperature and precipitation). We expect variability to increase with increased range size. To correct for this association, and to be able to include both variables in the model, we took the residuals from GAM models of the relationship between PCA components and range size.

Table S1: Loadings for the first two PCs, from a PCA on the variation (standard error, where n ~ 1,000) of each bioclim variable except the standard deviation of temperature seasonality (sd.bioclim4) and the standard error of precipitation seasonality (sd.bioclim15), as these are already measures of variability and based on other bioclimatic variables.

restricted.data <- complete.dataframe %>% drop_na(bioclim1) #Drop species without environmental data
PCA.bioclim <- prcomp(restricted.data[c('sd.bioclim1', 'sd.bioclim2', 'sd.bioclim3', 'sd.bioclim5', 'sd.bioclim6', 'sd.bioclim7', 'sd.bioclim8', 'sd.bioclim9', 'sd.bioclim10', 'sd.bioclim11', 'sd.bioclim12', 'sd.bioclim13', 'sd.bioclim14', 'sd.bioclim16', 'sd.bioclim17', 'sd.bioclim18', 'sd.bioclim19')], #Select sd.bioclim
                      scale = TRUE, center = TRUE)
PCA.predictions <- predict(PCA.bioclim)
restricted.data <- cbind(restricted.data, PCA.predictions)

#Create table with PC1 and PC2
as.data.frame(PCA.bioclim$rotation[,1:2]) %>% `colnames<-`(c("PC1", "PC2")) %>% pander()
  PC1 PC2
sd.bioclim1 -0.3323 0.08404
sd.bioclim2 -0.2386 0.05478
sd.bioclim3 -0.2323 -0.0647
sd.bioclim5 -0.3117 0.0818
sd.bioclim6 -0.3283 0.08438
sd.bioclim7 -0.2929 0.09912
sd.bioclim8 -0.2912 0.119
sd.bioclim9 -0.3181 0.1142
sd.bioclim10 -0.3108 0.06796
sd.bioclim11 -0.329 0.08941
sd.bioclim12 -0.1404 -0.4239
sd.bioclim13 -0.1608 -0.3072
sd.bioclim14 -0.01294 -0.3861
sd.bioclim16 -0.1721 -0.3242
sd.bioclim17 -0.02086 -0.3935
sd.bioclim18 -0.1448 -0.3057
sd.bioclim19 -0.01972 -0.3819
#run GAM of association between range size and spatial climatic variation
PC1.gam <- mgcv::gam(PC1*(-1) ~ s(log(range.size.m2)), data = restricted.data, family = "gaussian") #Inverse as the PC1 loads to the negative (COUNTER-INTUITIVE)

#take residuals
restricted.data$residuals.PC1 <- residuals.gam(PC1.gam) 

#We can do the same for PC2
PC2.gam <- mgcv::gam(PC2 ~ s(log(range.size.m2)), data = restricted.data, family = "gaussian")

#Take residuals
restricted.data$residuals.PC2 <- residuals.gam(PC2.gam)

#plot
par(mfrow = c(1,2))
plot.gam(PC1.gam, residuals = T, main = 'PC1 residuls (temp)')
plot.gam(PC2.gam, residuals = T, main = 'PC2 residuals (precip)')

Figure S1: The relationship between spatial variability in the temperature components (PC1) and log-range size is relatively strong. The relationship between spatial variability in precipitation (PC2) and log-range size is weaker. In the analysis we only used the residuals from PC1. PC1 accounts for 48.1 % of the variation.

Long term climate variability (LIG anomaly)

Climate stability through time can potentially affect diversification dynamics. To gain estimates for change in climate over the past ~130,000 years we used the difference in bioclim variables between the LIG and present values. The plots show the two PCs, broadly representing temperature and precipitation. We used the difference between the present and LIG climates as these represent a longer (evolutionarily meaningful) time-scale than the difference between the LGM and present climates.

Table S2: Loadings for the first two PCs of each of the PCAs for the difference in bioclimatic variables between today and the LIG. Here, PC1 is more heavily loaded for absolute temperature difference, while PC2 is more heavily loaded for absolute difference in precipitation.

historical.variation.data <- as.data.frame(restricted.data[1])

#FOR LIG

historical.variation.data$bio1.LIG.diff <- abs(restricted.data$bioclim1 - restricted.data$LIG.bi1)
historical.variation.data$bio2.LIG.diff <- abs(restricted.data$bioclim2 - restricted.data$LIG.bi2)
historical.variation.data$bio3.LIG.diff <- abs(restricted.data$bioclim3 - restricted.data$LIG.bi3)
historical.variation.data$bio4.LIG.diff <- abs(restricted.data$bioclim4 - restricted.data$LIG.bi4)
historical.variation.data$bio5.LIG.diff <- abs(restricted.data$bioclim5 - restricted.data$LIG.bi5)
historical.variation.data$bio6.LIG.diff <- abs(restricted.data$bioclim6 - restricted.data$LIG.bi6)
historical.variation.data$bio7.LIG.diff <- abs(restricted.data$bioclim7 - restricted.data$LIG.bi7)
historical.variation.data$bio8.LIG.diff <- abs(restricted.data$bioclim8 - restricted.data$LIG.bi8)
historical.variation.data$bio9.LIG.diff <- abs(restricted.data$bioclim9 - restricted.data$LIG.bi9)
historical.variation.data$bio10.LIG.diff <- abs(restricted.data$bioclim10 - restricted.data$LIG.bi10)
historical.variation.data$bio11.LIG.diff <- abs(restricted.data$bioclim11 - restricted.data$LIG.bi11)
historical.variation.data$bio12.LIG.diff <- abs(restricted.data$bioclim12 - restricted.data$LIG.bi12)
historical.variation.data$bio13.LIG.diff <- abs(restricted.data$bioclim13 - restricted.data$LIG.bi13)
historical.variation.data$bio14.LIG.diff <- abs(restricted.data$bioclim14 - restricted.data$LIG.bi14)
historical.variation.data$bio15.LIG.diff <- abs(restricted.data$bioclim15 - restricted.data$LIG.bi15)
historical.variation.data$bio16.LIG.diff <- abs(restricted.data$bioclim16 - restricted.data$LIG.bi16)
historical.variation.data$bio17.LIG.diff <- abs(restricted.data$bioclim17 - restricted.data$LIG.bi17)
historical.variation.data$bio18.LIG.diff <- abs(restricted.data$bioclim18 - restricted.data$LIG.bi18)
historical.variation.data$bio19.LIG.diff <- abs(restricted.data$bioclim19 - restricted.data$LIG.bi19)

historical.variation.data <- historical.variation.data %>% drop_na(bio1.LIG.diff)

#Run a PCA of the difference removing the variation variables (4 and 15)

#For LIG
PCA.LIG.bioclim <- prcomp(historical.variation.data[c(
'bio1.LIG.diff',
'bio2.LIG.diff',
'bio3.LIG.diff',
'bio5.LIG.diff',
'bio6.LIG.diff',
'bio7.LIG.diff',
'bio8.LIG.diff',
'bio9.LIG.diff',
'bio10.LIG.diff',
'bio11.LIG.diff',
'bio12.LIG.diff',
'bio13.LIG.diff',
'bio14.LIG.diff',
'bio16.LIG.diff',
'bio17.LIG.diff',
'bio18.LIG.diff',
'bio19.LIG.diff')],
                      scale = TRUE, center = TRUE)

#Create table with PC1 and PC2
as.data.frame(PCA.LIG.bioclim$rotation[,1:2]) %>% `colnames<-`(c("PC1.LIG", "PC2.LIG")) %>% pander()
  PC1.LIG PC2.LIG
bio1.LIG.diff -0.2956 -0.01943
bio2.LIG.diff -0.2442 0.0745
bio3.LIG.diff -0.1436 0.2085
bio5.LIG.diff -0.2867 -0.2566
bio6.LIG.diff -0.4094 0.03346
bio7.LIG.diff -0.4012 -0.09706
bio8.LIG.diff 0.07108 -0.2385
bio9.LIG.diff -0.3365 -0.00753
bio10.LIG.diff -0.07014 -0.3766
bio11.LIG.diff -0.4201 -0.004
bio12.LIG.diff -0.05181 0.3337
bio13.LIG.diff -0.2552 0.2591
bio14.LIG.diff 0.05701 0.3298
bio16.LIG.diff -0.1946 0.3164
bio17.LIG.diff 0.03896 0.3603
bio18.LIG.diff 0.06173 0.277
bio19.LIG.diff 0.08743 0.2856
#Now we can predict the PCA results: 
PCA.LIG.predictions <- as.data.frame(predict(PCA.LIG.bioclim)*(-1)) #So that higher numbers mean more variation
PCA.LIG.predictions <- rename(PCA.LIG.predictions, PC1.LIG = PC1,
       PC2.LIG = PC2)

#Bind them to dataframe, taking only the first two PCAs
historical.variation.data <- cbind(historical.variation.data, PCA.LIG.predictions[1:2])

#Now to the restricted dataframe
restricted.data <- right_join(restricted.data, historical.variation.data %>% dplyr::select(binomial, PC1.LIG, PC2.LIG), by = 'binomial')

PCA.plot <- historical.variation.data %>% ggplot(aes(x = PC1.LIG, y= PC2.LIG))+
  geom_point(shape = 21)+
  theme_minimal()

PCA.plot.m <- ggExtra::ggMarginal(
  p = PCA.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = 'gray'
)

grid.arrange(PCA.plot.m)

Figure S2: The distribution of the two PCA components for the absolute difference in bioclimatic variables between today and the LIG.

Correlations between environmental predictors

We checked the correlations between environmental predictors used in subsequent phylogenetic comparative analyses (PGLS models). Specifically we tested whether the following environmental variables are correlated:

  • Range size
  • Short temporal variability of temperature (mean BIOCLIM4)
  • Spatial variability of temperature (PCA1) [residual.PCA1]
  • Long-term variability of temperature (LIG)
  • NPP

To check the correlations between the environmental predictors to be used in the PGLS models we inspected the following correlation plots with the correlation value plotted:

#Draw plot
restricted.data$log.range.size <- log(restricted.data$range.size.m2)
pairs(restricted.data[,c("log.range.size", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP")], lower.panel=panel.smooth, upper.panel=panel.cor)

Figure S3: The highest correlation is between Long term temperature variation and seasonal variation. However, this correlation is moderate (0.53) so we can be confident that collinearity with not be an issue for our PGLS models.

Sexual Dimorphism

The sexual dimorphism dataset looks to have an over-dispersed distribution (see phylogenic plot in manuscript). Unfortunately, transformations do not greatly improve the distribution. Our model fit is expected to be reduced by this but further transformations are not obvious and would risk problems of interpretation. In any case, the phylogenetically corrected models we employ do not assume normality in the response variable.


Analysis

Simple speciation measures: Diversification rate (DR) and Node Density (ND)

We obtained two tip-rate metrics of speciation using statistics derived from the properties of the nodes and branches along root-to-tip paths of the phylogeny. Node density (ND) is a simple statistic calculating the density of nodes from the phylogenetic root-to-tip, while the diversification rate (DR) (e.g., Jetz et al. 2012; Quintero and Jetz 2018; Rabosky et al. 2018), is derived from the sum of edge lengths branching from a node, with each more basal node having the sum of lengths down-weighted.

The functions that generate these tip rates estimates are: calculate.log.es, calculate.nd and calculate.tb. Here they are run on 100 trees plus a full MCC tree of passerine birds based on 2500 sample trees of the posterior (Jetz et al. 2012). This was done using the phanghorn package (Schliep 2010).

#Obtain DR (log(es)) estimates, by calling the calculate.dr function
passerine.trees <- read.nexus('data/trees/Passerine_Trees_Full.nex')

#DR, ND and TB can be calculated on the 100 trees by: 
  
# es.list <- sapply(passerine.trees, calculate.log.es)
# nd.list <- sapply(passerine.trees, calculate.nd)
# tb.list <- sapply(passerine.trees, calculate.tb)

#To avoid re-running the time-consuming rates, we can save them and load them: 
es.list <- readRDS("data/es.list.rds")
nd.list <- readRDS("data/nd.list.rds")
tb.list <- readRDS("data/tb.list.rds")

#Also load in our MCC tree
MCC.passerine <- read.tree('data/MCC.passerine.tre') #Tree based on 2500 samples of posterior

To gain an estimate of variation between trees we can take the values and obtain a mean, variance and coefficient of variation (CV). We obtain a CV value to account for unequal variation at higher estimates. The CV is expressed as the ratio of the standard deviation (\(\sigma\)) to the mean (\(\mu\)) (converting to a percentage not needed):

\(C_V = \frac{\sigma}{\mu} \times 100\)

Furthermore, we can account for over-dispersion by obtaining a transformed version of CV. In such a case we can calculate CV for a log-normal distribution:

\({\widehat {c_{\rm {v}}}}_{\rm {ln}}={\sqrt {\mathrm {e} ^{s_{\rm {ln}}^{2}}-1}}\)

Given that our ES value is log-transformed, the calculation of \({\widehat {c_{\rm {v}}}}_{\rm {ln}}\) is perhaps more ideal. As for ND we can see from the a qqnorm plot that value of ND is normally distributed, therefore the regular CV is more appropriate for this distribution.

#For ND
#transpose the list so that the elements are the species
t.nd.list <- nd.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
nd.values <- as.data.frame(t(sapply(t.nd.list, function(x) sapply(x, function(x) (x)))))
nd.values <- nd.values %>% tibble::rownames_to_column("binomial")
nd.summary <- melt(nd.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.nd = mean(value), #use normal mean, not a rate
  var.nd = var(value), #Don't really need this if we use CV
  CV.nd = (sd(value)/mean(value)) #Given normal distribution of ND, standard CV formula is appropriate
)

#For ES
#transpose the list so that the elements are the species
t.es.list <- es.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
es.values <- as.data.frame(t(sapply(t.es.list, function(x) sapply(x, function(x) (x)))))
es.values <- es.values %>% tibble::rownames_to_column("binomial")
es.summary <- melt(es.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.loges = log(1/mean(1/exp(value))), #use harmonic mean for rates and log-transform post calculation of HM
  var.loges = var(value), #Don't really need this if we use CV
  CV.loges = sqrt(exp(var(value)) - 1) #CV for log-normal distribution of ES
)

#For TB
#transpose the list so that the elements are the species
t.tb.list <- tb.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
tb.values <- as.data.frame(t(sapply(t.tb.list, function(x) sapply(x, function(x) (x)))))
tb.values <- tb.values %>% tibble::rownames_to_column("binomial")
tb.summary <- melt(tb.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.tb = mean(value), #use normal mean; not a rate
  var.tb = var(value), #Don't really need this if we use CV
  CV.tb = (sd(value)/mean(value)) #Given normal distribution of tb, standard CV formula is appropriate
)

MCC.nd.df <- as.data.frame(calculate.nd(MCC.passerine)) %>% tibble::rownames_to_column("binomial")
MCC.dr.df <- as.data.frame(calculate.log.es(MCC.passerine)) %>% tibble::rownames_to_column("binomial")

#Join the dataframes
dr.summary <- full_join(nd.summary, es.summary, by = "binomial")
dr.summary <- full_join(dr.summary, tb.summary, by = "binomial")
dr.summary <- full_join(dr.summary, MCC.nd.df, by = "binomial")
dr.summary <- full_join(dr.summary, MCC.dr.df, by = "binomial")

dr.summary <- dr.summary %>% rename(TipLabel = binomial, MCC.ND = nd, MCC.DR = es)

saveRDS(dr.summary, 'data/dr.summary.rds')
dr.summary <- readRDS('data/dr.summary.rds')
#Plot a summary of logES and TB with the DRs with their weightings
dr.plot <- dr.summary %>% ggplot(aes(x = mean.loges, y = mean.nd, 
                          fill = (1/(CV.nd))/sum((1/(CV.nd)))*100, 
                      size = (1/(CV.loges))/sum((1/(CV.loges)))*100))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = mean.loges - 0.674*sqrt(var.loges), xmax = mean.loges + 0.674*sqrt(var.loges)), 
                 size = 0.015)+
  geom_errorbar(aes(ymin = mean.nd - 0.674*sqrt(var.nd), ymax = mean.nd + 0.674*sqrt(var.nd)), 
                size = 0.015)+
  theme(legend.position="bottom")+
  labs(size = 'Weight (%) [using DR]', y='Node Density [ND]', x= 'Log Diversification Rate [DR]', fill = 'Weight (%) [using ND] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "RdPu", direction = 1, breaks=c(0,0.015,.03),
                           limits=c(0,0.03))
marginal.dr <-ggMarginal(
  p = dr.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = '#0000009C',
  fill = '#D12E6769'
)
grid.newpage()
grid.draw(marginal.dr)

Figure S4: \(\lambda_{DR}\) and \(\lambda_{ND}\) estimates of speciation rate are correlated. Furthermore, we see similar variation across each species in the estimates across 100 trees as the inverse of the coefficient of variation (weight) is approximately equivalent across all values of ND and DR (colours and size). Notably, these estimates speciation rate are relatively uncertain across all species, with 50 % CIs in both axes plotted.

restricted.data <- left_join(restricted.data, dr.summary, by = 'TipLabel')

es.plot <- restricted.data %>% ggplot(aes(y = mean.loges, x = SDi, fill = mean.loges))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "YlOrBr", direction = 1, guide = FALSE)+
ylab("mean.DR")+
theme_minimal()

nd.plot <- restricted.data %>% ggplot(aes(y = mean.nd, x = SDi, fill = mean.nd))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "Greens", direction = 1, guide = FALSE)+
theme_minimal()

tb.plot <- restricted.data %>% ggplot(aes(y = mean.tb, x = SDi, fill = mean.tb))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "PuBu", direction = 1, guide = FALSE)+
theme_minimal()

grid.arrange(es.plot, nd.plot, tb.plot, nrow = 3)

Figure S5: Scatter plots showing the raw relationship between sexual dichromatism (SDi) and speciation rates. Across all three measures of speciation the pattern and spread is similar, with no obvious relationship. Note that terminal branch length (mean.tb) was not used in the analysis.

BAMM measures of speciation and extinction

Set up BAMM parameters

For the use of BAMM we used the following code to generate the parameters across the 100 trees. Each parameter value is specified in this code chunk to ensure reproducibility. These same parameters were also used for the MCC tree (with the seed set at 2500).

name.passerine.tree <- names(passerine.trees)

priors <- sapply(name.passerine.tree, function(x) {
  setBAMMpriors(passerine.trees[[x]], outfile = NULL)
})

sapply(name.passerine.tree, function(x) {
  write.tree(passerine.trees[[x]], paste("data/bamm_files/", x, ".tre", sep=""))
})

# Here is a block of parameters for the control file. We can make a control file for each tree:
params <- list()
for (x in name.passerine.tree) {

# GENERAL SETUP AND DATA INPUT

params[[x]] <- list(modeltype = 'speciationextinction',
# Specify "speciationextinction" or "trait" analysis
                                  
treefile = paste(x, ".tre", sep=""),
# File name of the phylogenetic tree to be analyzed

runInfoFilename = 'run_info.txt',
# File name to output general information about this run

sampleFromPriorOnly = 0,
# Whether to perform analysis sampling from prior only (no likelihoods computed)

runMCMC = 1,
# Whether to perform the MCMC simulation. If runMCMC = 0, the program will only
# check whether the data file can be read and the initial likelihood computed

loadEventData = 0,                       
# Whether to load a previous event data file

eventDataInfile = 'event_data_in.txt',
# File name of the event data file to load, used only if loadEventData = 1

initializeModel = 1,
# Whether to initialize (but not run) the MCMC. If initializeModel = 0, the
# program will only ensure that the data files (e.g., treefile) can be read

useGlobalSamplingProbability = 1,
# Whether to use a "global" sampling probability. If False (0), expects a file
# name for species-specific sampling probabilities (see sampleProbsFilename)
                                        
globalSamplingFraction = 1,
# The sampling probability. If useGlobalSamplingProbability = 0, this is ignored
# and BAMM looks for a file name with species-specific sampling fractions

sampleProbsFilename = 'sample_probs.txt',
# File name containing species-specific sampling fractions

seed = as.numeric(gsub("tree_", "", x, perl = TRUE)),
# Seed for the random number generator. Set for reproducibility to the number of the treefile

overwrite = 1,
# If True (1), the program will overwrite any output files in the current
# directory (if present)


# PRIORS

expectedNumberOfShifts = 100,
# prior on the number of shifts in diversification
# Suggested values: 
#     expectedNumberOfShifts = 1.0 for small trees (< 500 tips)
#  expectedNumberOfShifts = 10 or even 50 for large trees (> 5000 tips) 
 
lambdaInitPrior = as.numeric(priors['lambdaInitPrior', x]),
# Prior (rate parameter of exponential) on the initial lambda value for rate
# regimes

lambdaShiftPrior = 0.05,
# Prior (std dev of normal) on lambda shift parameter for rate regimes
# You cannot adjust the mean of this distribution (fixed at zero, which is
# equal to a constant rate diversification process)

muInitPrior = as.numeric(priors['muInitPrior', x]),
# Prior (rate parameter of exponential) on extinction rates  

lambdaIsTimeVariablePrior = 1,
# Prior (probability) of the time mode being time-variable (vs. time-constant)
            

# MCMC SIMULATION SETTINGS & OUTPUT OPTIONS

numberOfGenerations = '100000000',
# Number of generations to perform MCMC simulation

mcmcOutfile = 'mcmc_out.txt',
# File name for the MCMC output, which only includes summary information about
# MCMC simulation (e.g., log-likelihoods, log-prior, number of processes)

mcmcWriteFreq = 1000,
# Frequency in which to write the MCMC output to a file

eventDataOutfile = 'event_data.txt',
# The raw event data (these are the main results). ALL of the results are
# contained in this file, and all branch-specific speciation rates, shift
# positions, marginal distributions etc can be reconstructed from this output.
# See R package BAMMtools for working with this output

eventDataWriteFreq = 1000,
# Frequency in which to write the event data to a file

printFreq = 10000,
# Frequency in which to print MCMC status to the screen

acceptanceResetFreq = 1000,
# Frequency in which to reset the acceptance rate calculation
# The acceptance rate is output to both the MCMC data file and the screen

outName = x,
# Optional name that will be prefixed on all output files (separated with "_")
# If commented out, no prefix will be used


# OPERATORS: MCMC SCALING OPERATORS

updateLambdaInitScale = 2,
# Scale parameter for updating the initial speciation rate for each process

updateLambdaShiftScale = 0.1,
# Scale parameter for the exponential change parameter for speciation

updateMuInitScale = 2,
# Scale parameter for updating initial extinction rate for each process

updateEventLocationScale = 0.1,
# Scale parameter for updating LOCAL moves of events on the tree
# This defines the width of the sliding window proposal
 
updateEventRateScale = 4,
# Scale parameter (proportional shrinking/expanding) for updating
# the rate parameter of the Poisson process

# OPERATORS: MCMC MOVE FREQUENCIES

updateRateEventNumber = 1,
# Relative frequency of MCMC moves that change the number of events

updateRateEventPosition = 0.25,
# Relative frequency of MCMC moves that change the location of an event on the
# tree

updateRateEventRate = 1,
# Relative frequency of MCMC moves that change the rate at which events occur 

updateRateLambda0 = 1,
# Relative frequency of MCMC moves that change the initial speciation rate
# associated with an event

updateRateLambdaShift = 1,
# Relative frequency of MCMC moves that change the exponential shift parameter
# of the speciation rate associated with an event

updateRateMu0 = 1,
# Relative frequency of MCMC moves that change the extinction rate for a given
# event

updateRateLambdaTimeMode = 0,
# Relative frequency of MCMC moves that flip the time mode
# (time-constant <=> time-variable)

localGlobalMoveRatio = 10,
# Ratio of local to global moves of events 


# INITIAL PARAMETER VALUES

lambdaInit0 = 0.032,
# Initial speciation rate (at the root of the tree)

lambdaShift0 = 0,
# Initial shift parameter for the root process

muInit0 = 0.005,
# Initial value of extinction (at the root)

initialNumberEvents = 0,
# Initial number of non-root processes


# METROPOLIS COUPLED MCMC

numberOfChains = 1,
# Number of Markov chains to run

deltaT = 0.01,
# Temperature increment parameter. This value should be > 0
# The temperature for the i-th chain is computed as 1 / [1 + deltaT * (i - 1)]

swapPeriod = 1000,
# Number of generations in which to propose a chain swap

chainSwapFileName = 'chain_swap.txt',
# File name in which to output data about each chain swap proposal.
# The format of each line is [generation],[rank_1],[rank_2],[swap_accepted]
# where [generation] is the generation in which the swap proposal was made,
# [rank_1] and [rank_2] are the chains that were chosen, and [swap_accepted] is
# whether the swap was made. The cold chain has a rank of 1.


# NUMERICAL AND OTHER PARAMETERS

minCladeSizeForShift = 3,
# Allows you to constrain location of possible rate-change events to occur
# only on branches with at least this many descendant tips. A value of 1
# allows shifts to occur on all branches. 

segLength = 0.025,
# Controls the "grain" of the likelihood calculations. Approximates the
# continuous-time change in diversification rates by breaking each branch into
# a constant-rate diversification segments, with each segment given a length
# determined by segLength. segLength is in units of the root-to-tip distance of
# the tree. So, if the segLength parameter is 0.01, and the crown age of your
# tree is 50, the "step size" of the constant rate approximation will be 0.5.
# If the value is greater than the branch length (e.g., you have a branch of
# length < 0.5 in the preceding example) BAMM will not break the branch into
# segments but use the mean rate across the entire branch.

outName = x)
  }

bammcontrolfile <- list()
for (x in name.passerine.tree) {
  bammcontrolfile[x] <- paste("data/bamm_files/control_", x, ".txt", sep="")
}

# Now writing control parameters to file

for (x in name.passerine.tree) {generateControlFile(file = bammcontrolfile[[x]], type = "diversification", params = params[[x]])}

Run analysis

BAMM can be run through the terminal through the following syntax: bamm -c control_tree_xxxx.txt. To generate these commands we can use a loop function, from which we get:

bamm.commands <- list()
for (x in name.passerine.tree) {
  bamm.commands[x] <- paste("bamm -c control_", x, ".txt", sep="")
}

The analysis was run over multiple CPU’s, each generating a respective MCMC and EventData output. Due to the size of the event data file (~ 50 Gb in total) they are not included as supplementary material here. However we have simplified the event data objects into tip rate estimates of the mean and variance across 100 trees + MCC.

Read in the event data and extract tip data

Now we have a series of trees and event data we can read in. Firstly and most importantly let’s check for convergence in the MCC tree through the following:

Table S3: Effective sample size (ESS) for the two key BAMM parameters (number of evolutionary shifts and log-Likelihood) for the run on the MCC indicate that BAMM converges (ESS > 200).

#Read in the tree and MCMC to check for convergence 
MCC.BAMM.tree  <- read.tree("data/BAMM_MCC/MCC.passerine.tre") #Same as other MCC tree already loaded
MCC.BAMM.mcmc <- read.csv( "data/BAMM_MCC/tree_MCC_mcmc_out.txt" , stringsAsFactors=F)

#Plot of convergence can be generated by:
#plot(MCC.BAMM.mcmc$logLik ~ MCC.BAMM.mcmc$generation)

#Looks like it has converged so let's discard burn in: 
burnstart <- floor(0.1 * nrow(MCC.BAMM.mcmc))
postburn <- MCC.BAMM.mcmc[burnstart:nrow(MCC.BAMM.mcmc), ]

#We can also check effective population sizes of the log-likelihhod and number of shift events in each sample
#We want at least 200 (although that's on the low side)

cbind(effectiveSize(postburn$N_shifts), effectiveSize(postburn$logLik)) %>% `colnames<-`(c("N_Shifts", "logLik")) %>% `rownames<-`("Effective Sample Size") %>% pander()
  N_Shifts logLik
Effective Sample Size 2456 1302

Those values are well above 200, so we can use the MCC data in analysis! We can also check the convergence of BAMM across 100 runs of BAMM. The Raw MCMCs are not included in this file or repository but we can read in a dataframe that has extracted the ESS for each of the runs.

Table S4: For the 100 trees that BAMM was run on effective sample size (ESS) for the two key BAMM parameters (number of evolutionary shifts and log-Likelihood) also indicates that BAMM converges with the minimum for each parameter across the 100 trees being over 200.

ESS <- readRDS('data/ESS.rds')
summary(ESS) %>% pander()
N_Shifts logLik
Min. : 723.4 Min. : 278.7
1st Qu.:1396.1 1st Qu.: 511.3
Median :1618.6 Median : 638.6
Mean :1709.5 Mean : 669.9
3rd Qu.:2070.8 3rd Qu.: 804.8
Max. :3425.6 Max. :1244.1

Given that it appears BAMM converges we can we can make a data frame with the mean and variance for extinction and speciation tip rates from the large event data set for the MCC with the following code and then plot the variation we see in tip-rates.

# Read in Event Data
MCC.BAMM.ED <- getEventData(MCC.BAMM.tree,  "data/BAMM_MCC/tree_MCC_event_data.txt", burnin=0.1, nsamples=1000)
## Reading event datafile:  data/BAMM_MCC/tree_MCC_event_data.txt 
##      ...........
## Read a total of 100000 samples from posterior
## 
## Discarded as burnin: GENERATIONS <  9999000
## Analyzing  1000  samples from posterior
## 
## Setting recursive sequence on tree...
## 
## Done with recursive sequence
saveRDS(MCC.BAMM.ED, 'data/MCC.BAMM.ED.rds')
#From the Event Data we can extract 

library(purrr)
#BAMM.EventData <- readRDS('data/BAMM.EventData.rds')

#Big lapply over each tree in BAMM event data
BAMM.extraction.function <- function(x) {
######Get mean and var for lambda
#Transpose list so each element in the list is a species
transposed.lambda <- lapply(purrr::transpose(x$tipLambda), unlist)

#Now turn it into a df with mean and variance
lambda <- sapply(transposed.lambda, function(x) {
  mean.lambda = mean(log(x))
  var.lambda = var(log(x))
  return(c(mean.lambda, var.lambda))
}) %>% t() %>% as.data.frame() %>% `colnames<-`(c("mean.lambda", "var.lambda"))

lambda$TipLabel <- x[["tip.label"]]

#####NOW FOR Extinction

#Transpose list so each element in the list is a species
transposed.mu <- lapply(purrr::transpose(x$tipMu), unlist)

#Now turn it into a df with mean and variance
mu <- sapply(transposed.mu, function(x) {
  mean.mu = mean(log(x))
  var.mu = var(log(x))
  return(c(mean.mu, var.mu))
}) %>% t() %>% as.data.frame() %>% `colnames<-`(c("mean.mu", "var.mu"))

mu$TipLabel <- x[["tip.label"]]

left_join(lambda, mu, by = "TipLabel")
}

MCC.BAMM.df <- BAMM.extraction.function(MCC.BAMM.ED)

#Save df for later use
saveRDS(MCC.BAMM.df, 'data/MCC.BAMM.df.rds')
MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')
#Plot a summary of logES and TB with the DRs with their weightings
BAMM.MCC.plot <- MCC.BAMM.df %>% ggplot(aes(x = mean.lambda, y = mean.mu, 
                          fill = 1/var.lambda, 
                      size = 1/var.mu))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = mean.lambda - 0.674*sqrt(var.lambda), xmax = mean.lambda + 0.674*sqrt(var.lambda)),
                 size = 0.0025)+
  geom_errorbar(aes(ymin = mean.mu - 0.674*sqrt(var.mu), ymax = mean.mu + 0.674*sqrt(var.mu)),
                size = 0.0025)+
  # scale_y_continuous(trans = "log")+
  # scale_x_continuous(trans = "log")+
  # xlim(-8, 3)+
  # ylim(-10,2)+
  # scale_x_continuous(trans = "log")+
  # scale_y_continuous(trans = "log")+
  theme(legend.position="bottom")+
  labs(size = 'Inverse log(var) / Weight [using Lambda]', y='Log Extinction [Mu]', x= 'Log Speciation [Lambda]', fill = 'Inverse log(var) / Weight [using Mu] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "Reds", direction = 1)

BAMM.variance <- ggExtra::ggMarginal(
  p = BAMM.MCC.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = '#BA3B1C91'
)
grid.newpage()
grid.draw(BAMM.variance)

Figure S6: The tip-rate estimates for BAMM are highly variable within each run of BAMM. Across most species there is high variability in the posterior distribution of tip-rate estimates. Here we show mean values, weights based on the variance and 50 % CIs.

Analysis of BAMM results

Given the high variability in tip-rate estimates from the above plot, below we peformed some diagnostics on BAMM to demonstrate that the variability is unlikely an error in sampling or parameters, rather an inherent aspect of BAMM. The following is not inherently necessary to understand the conclusions drawn from the paper, however they do raise a set of methodological questions about BAMM that warrant further investigation.

Credible number of shifts

To plot the credible shift set, we need the prior distribution on the number of rate shifts (this is generated internally by BAMMtools). We can then estimate the credible set of rate shifts using the BAMMtools function credibleShiftSet:

css <- credibleShiftSet(MCC.BAMM.ED, expectedNumberOfShifts=100, threshold=5, set.limit = 0.95)

#Now we obtain the number of distinct shifts: (Out of 1000 samples this is super high, essentially each one is distinct)

summary(css)
## 
##  95 % credible set of rate shift configurations sampled with BAMM
## 
## Distinct shift configurations in credible set:  950
## 
## Frequency of 9 shift configurations with highest posterior probability:
## 
## 
##    rank     probability cumulative  Core_shifts
##          1      0.001      0.001         54
##          2      0.001      0.002         44
##          3      0.001      0.003         46
##          4      0.001      0.004         47
##          5      0.001      0.005         49
##          6      0.001      0.006         45
##          7      0.001      0.007         47
##          8      0.001      0.008         43
##          9      0.001      0.009         45
## 
## ...omitted 941 additional distinct shift configurations
## from the credible set. You can access the full set from your 
## credibleshiftset object

Notably, each of the top shifts has low probability, indicating that out of the entire posterior sample we cannot differentiate which shift configuration is more likely. Based on the following estimates of the BayesFactor, we are confident that the number of shifts is non-zero, however it is still quite a wide distribution.

round(computeBayesFactors(MCC.BAMM.mcmc, expectedNumberOfShifts=100, burnin=0.1)[,1], digits = 2)
##      42      43      44      45      46      47      48      49      50 
##    1.00   17.17   24.48   49.45  120.71  283.77  454.33  784.80 1291.85 
##      51      52      53      54      55      56      57      58      59 
## 1813.33 2561.62 3608.07 4574.91 5492.44 6322.11 7141.12 7720.26 7886.28 
##      60      61      62      63      64      65      66      67      68 
## 7846.73 7661.83 7116.15 6466.36 5749.34 4811.16 4136.80 3338.17 2674.70 
##      69      70      71      72      73      74      75      76      77 
## 2137.61 1593.48 1218.40  868.01  643.91  452.36  311.07  199.17  135.99 
##      78      79      80      81      82      83      84      85      87 
##   85.85   65.03   40.87   14.74   13.40    9.02   10.63    1.53    1.56
plotPrior(MCC.BAMM.mcmc, expectedNumberOfShifts=100)

Figure S7: The apparent convergence of the number of shifts in the posterior is at odds with the variability seen in the CSS and although there is greater certainty in the number of shifts being within the range above, the position of those shifts remains variable.

We can compare our run of BAMM against Harvey et al. (2017) who used BAMM on a genetic-only MCC tree with different parameters.

load('data/Hackett_split_eventsample.rda')
css3 <- credibleShiftSet(ed, expectedNumberOfShifts=100, threshold=5, set.limit = 0.95)

Harvey.BAMM.df <- BAMM.extraction.function(ed)

Harvey.BAMM.df %>% ggplot(aes(x = mean.lambda, y = mean.mu, 
                          fill = 1/(var.lambda), 
                      size = 1/(var.mu)))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = mean.lambda - 0.674*sqrt(var.lambda), xmax = mean.lambda + 0.674*sqrt(var.lambda)), 
                 size = 0.0025)+
  geom_errorbar(aes(ymin = mean.mu - 0.674*sqrt(var.mu), ymax = mean.mu + 0.674*sqrt(var.mu)), 
                size = 0.0025)+
  # scale_x_continuous(trans = "log")+
  # scale_y_continuous(trans = "log")+
  theme(legend.position="bottom")+
  labs(size = 'Weight [using Lambda]', y='Log Extinction [Mu]', x= 'Log Speciation [Lambda]', fill = 'Weight [using Mu] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "Reds", direction = 1)

# BAMM.variance <- ggExtra::ggMarginal(
#   p = BAMM.MCC.plot,
#   type = 'density',
#   margins = 'both',
#   size = 5,
#   colour = 'black',
#   fill = '#BA3B1C91'
# )

Figure S8: Based on the event data from the BAMM run by Harvey et al. (2017) we find that in both our case and theirs the tip-rate estimates for many species are extremely variable across samples of the posterior probability. Here we present mean estimates with 50 % CIs.

PGLS Models

The method behind running the PGLS models is as follows:

1: Estimate the phylogenetic signal (\(\lambda\)) by running a model without interactions and all six predictor variables. The value of \(\lambda\) obtained here will be fixed in all successive models. The value is fixed for subsequent models as independently estimating it in each case becomes computationally intensive.

2: Create a global model of six predictor variables plus the five interactions between sexual dichromatism and the other variables.

3: Dredge the global model but fix the six independent predictors, hence conducting model selection on the interaction terms.

4: Take the top model in the MCC model and run it on the 100 phylogenetic trees.

5: Repeat this step for DR, ND, BAMM-speciation and BAMM-extinction.

PGLS Models on DR and ND

Using corPagel we can estimate the phylogenetic signal for a model with all predictors (interactions do not appear to affect the estimate of \(\lambda\)):

#Prune tree
pruned.MCC.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(restricted.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(restricted.data) <- restricted.data$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.corPagel <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.tree, fixed = FALSE), 
                data = restricted.data, 
                method = "REML")
saveRDS(MCC.DR.corPagel, 'data/MCC.DR.corPagel.rds')

#Run a corPagel model to estimate lambda for ND
MCC.ND.corPagel <- gls(MCC.ND ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.tree, fixed = FALSE), 
                data = restricted.data, 
                method = "REML")
saveRDS(MCC.ND.corPagel, 'data/MCC.ND.corPagel.rds')

Inspect the \(\lambda\) value, which can then be fixed for successive models.

MCC.DR.corPagel <- readRDS('data/MCC.DR.corPagel.rds')
MCC.ND.corPagel <- readRDS('data/MCC.ND.corPagel.rds')
MCC.DR.corPagel[["modelStruct"]][["corStruct"]] %>% `names<-`("DR lambda") %>% pander()
  • DR lambda: 0.985
MCC.ND.corPagel[["modelStruct"]][["corStruct"]] %>% `names<-`("ND lambda") %>% pander()
  • ND lambda: 0.9996

The lambda is high and similar to if we assume Brownian Motion \(\lambda = 1\).However, given the large sample size, this difference may have an effect on the results so we included it as a fixed value for \(\lambda\) in all successive models.

#Run model for DR
MCC.DR <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                correlation = corPagel(0.985, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

#Run model for ND
MCC.ND <- gls(MCC.ND ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP
                      + SDi*log(range.size.m2)
                      + SDi*bioclim4
                      + SDi*residuals.PC1
                      + SDi*PC1.LIG
                      + SDi*NPP,
                correlation = corPagel(0.9996, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("restricted.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

dredged.ND.model <- pdredge(MCC.ND, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust, trace = TRUE)

saveRDS(dredged.ND.model, "data/dredged.ND.model.rds")

dredged.DR.model <- pdredge(MCC.DR, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust, trace = TRUE)

saveRDS(dredged.DR.model, "data/dredged.DR.model.rds")

Table S5: The dredged models both show the top model is one with no interactions, with \(\delta AICc > 4\) in both cases. We can be reasonably confident that interactions are unlikely to affect the pattern of speciation we see in passerine birds.

dredged.DR.model <- readRDS("data/dredged.DR.model.rds")
dredged.ND.model <- readRDS("data/dredged.ND.model.rds")
kable(dredged.DR.model, "html", caption = "DR Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
DR Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 SDi bioclim4:SDi log(range.size.m2):SDi NPP:SDi PC1.LIG:SDi residuals.PC1:SDi df logLik AICc delta weight
0 -2.879726 2.14e-05 -0.0065788 3e-07 0.0015918 -0.0024409 -0.0012788 NA NA NA NA NA 8 -3409.543 6835.112 0.00000 9.988662e-01
8 -2.879409 2.14e-05 -0.0065797 3e-07 -0.0000977 -0.0023890 -0.0012943 NA NA NA 0.0003488 NA 9 -3416.085 6850.201 15.08962 5.282489e-04
16 -2.879275 2.15e-05 -0.0066037 3e-07 0.0015571 -0.0016227 -0.0013077 NA NA NA NA -0.0001614 9 -3416.488 6851.007 15.89504 3.531362e-04
2 -2.884255 2.13e-05 -0.0064061 3e-07 0.0016045 -0.0024599 -0.0003044 NA -0.0000369 NA NA NA 9 -3416.839 6851.708 16.59642 2.486802e-04
1 -2.878346 1.57e-05 -0.0065697 3e-07 0.0015926 -0.0024296 -0.0015877 1.2e-06 NA NA NA NA 9 -3421.207 6860.445 25.33357 3.150599e-06
24 -2.878637 2.16e-05 -0.0066194 3e-07 -0.0004578 -0.0010799 -0.0013431 NA NA NA 0.0004118 -0.0002563 10 -3422.849 6865.736 30.62451 2.236036e-07
4 -2.884726 2.14e-05 -0.0065588 7e-07 0.0015614 -0.0024120 -0.0003592 NA NA -1e-07 NA NA 9 -3424.304 6866.638 31.52657 1.424293e-07
10 -2.890366 2.13e-05 -0.0061612 3e-07 -0.0001745 -0.0024316 0.0010660 NA -0.0000894 NA 0.0003710 NA 10 -3423.318 6866.674 31.56272 1.398780e-07
18 -2.885300 2.14e-05 -0.0063743 3e-07 0.0015727 -0.0016177 -0.0000092 NA -0.0000492 NA NA -0.0001673 10 -3423.771 6867.580 32.46864 8.892676e-08
9 -2.879392 2.13e-05 -0.0065796 3e-07 -0.0000944 -0.0023889 -0.0012984 0.0e+00 NA NA 0.0003481 NA 10 -3427.742 6875.522 40.41027 1.676986e-09
17 -2.876923 1.26e-05 -0.0065995 3e-07 0.0015442 -0.0012710 -0.0018051 1.9e-06 NA NA NA -0.0002272 10 -3428.011 6876.060 40.94823 1.281480e-09
3 -2.889344 1.34e-05 -0.0061267 3e-07 0.0016253 -0.0024734 0.0007735 1.7e-06 -0.0000939 NA NA NA 10 -3428.369 6876.775 41.66345 8.962004e-10
12 -2.884141 2.14e-05 -0.0065608 6e-07 -0.0000865 -0.0023628 -0.0004253 NA NA -1e-07 0.0003406 NA 10 -3430.868 6881.773 46.66161 7.363229e-11
26 -2.893368 2.14e-05 -0.0060577 3e-07 -0.0005938 -0.0010207 0.0018405 NA -0.0001207 NA 0.0004474 -0.0002792 11 -3430.029 6882.104 46.99264 6.240023e-11
20 -2.884058 2.15e-05 -0.0065821 7e-07 0.0015321 -0.0016854 -0.0004342 NA NA -1e-07 NA -0.0001436 10 -3431.268 6882.574 47.46202 4.934687e-11
6 -2.894916 2.13e-05 -0.0061873 7e-07 0.0015856 -0.0024497 0.0018171 NA -0.0000789 -1e-07 NA NA 10 -3431.542 6883.123 48.01090 3.750351e-11
25 -2.877781 1.82e-05 -0.0066168 3e-07 -0.0003387 -0.0009795 -0.0015303 7.0e-07 NA NA 0.0003864 -0.0002755 11 -3434.451 6890.947 55.83508 7.500203e-13
11 -2.891638 1.89e-05 -0.0060871 3e-07 -0.0000838 -0.0024371 0.0013317 5.0e-07 -0.0001044 NA 0.0003534 NA 11 -3434.891 6891.827 56.71540 4.829629e-13
5 -2.883337 1.59e-05 -0.0065503 7e-07 0.0015625 -0.0024015 -0.0006676 1.1e-06 NA -1e-07 NA NA 10 -3435.973 6891.983 56.87160 4.466790e-13
19 -2.894530 8.20e-06 -0.0058842 3e-07 0.0015867 -0.0010875 0.0019942 2.8e-06 -0.0001529 NA NA -0.0002772 11 -3435.062 6892.169 57.05691 4.071512e-13
28 -2.882932 2.16e-05 -0.0065996 6e-07 -0.0004216 -0.0011518 -0.0005610 NA NA -1e-07 0.0003998 -0.0002377 11 -3437.666 6897.377 62.26522 3.011531e-14
14 -2.901098 2.13e-05 -0.0059406 7e-07 -0.0001980 -0.0024212 0.0032021 NA -0.0001318 -1e-07 0.0003719 NA 11 -3438.019 6898.084 62.97271 2.114250e-14
22 -2.895431 2.14e-05 -0.0061672 7e-07 0.0015574 -0.0016836 0.0019991 NA -0.0000884 -1e-07 NA -0.0001523 11 -3438.493 6899.032 63.92041 1.316333e-14
27 -2.898245 1.33e-05 -0.0057878 3e-07 -0.0003739 -0.0007535 0.0028896 1.7e-06 -0.0001774 NA 0.0004038 -0.0003357 12 -3441.456 6906.965 71.85360 2.492851e-16
13 -2.884149 2.14e-05 -0.0065608 6e-07 -0.0000879 -0.0023629 -0.0004235 0.0e+00 NA -1e-07 0.0003409 NA 11 -3442.524 6907.094 71.98274 2.336973e-16
21 -2.881617 1.31e-05 -0.0065791 6e-07 0.0015211 -0.0013509 -0.0009438 1.8e-06 NA -1e-07 NA -0.0002065 11 -3442.807 6907.660 72.54812 1.761499e-16
7 -2.901185 1.25e-05 -0.0058634 7e-07 0.0016077 -0.0024642 0.0031373 1.8e-06 -0.0001448 -1e-07 NA NA 11 -3443.043 6908.132 73.02027 1.391094e-16
30 -2.903189 2.14e-05 -0.0058581 7e-07 -0.0005923 -0.0010894 0.0037855 NA -0.0001584 -1e-07 0.0004440 -0.0002637 12 -3444.762 6913.578 78.46668 9.134510e-18
29 -2.882110 1.85e-05 -0.0065975 6e-07 -0.0003149 -0.0010605 -0.0007389 6.0e-07 NA -1e-07 0.0003772 -0.0002552 12 -3449.270 6922.594 87.48267 1.006670e-19
15 -2.903111 1.78e-05 -0.0058318 7e-07 -0.0000710 -0.0024288 0.0036204 7.0e-07 -0.0001537 -1e-07 0.0003472 NA 12 -3449.580 6923.214 88.10213 7.385375e-20
23 -2.905738 7.50e-06 -0.0056398 7e-07 0.0015712 -0.0011302 0.0042249 2.9e-06 -0.0001998 -1e-07 NA -0.0002669 12 -3449.755 6923.565 88.45298 6.197088e-20
31 -2.908923 1.25e-05 -0.0055552 7e-07 -0.0003523 -0.0008007 0.0050118 1.9e-06 -0.0002219 -1e-07 0.0003963 -0.0003247 13 -3456.170 6938.402 103.29021 3.718126e-23
kable(dredged.ND.model, "html", caption = "ND Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
ND Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 SDi bioclim4:SDi log(range.size.m2):SDi NPP:SDi PC1.LIG:SDi residuals.PC1:SDi df logLik AICc delta weight
0 0.1404213 4e-07 -0.0001462 0 -0.0000102 -6.40e-05 -0.0000575 NA NA NA NA NA 8 14329.63 -28643.24 0.00000 9.999550e-01
8 0.1403572 4e-07 -0.0001421 0 -0.0000850 -5.65e-05 -0.0000618 NA NA NA 1.49e-05 NA 9 14320.04 -28622.05 21.19489 2.497865e-05
16 0.1404289 5e-07 -0.0001465 0 -0.0000108 -5.56e-05 -0.0000581 NA NA NA NA -1.7e-06 9 14319.24 -28620.45 22.79627 1.121590e-05
2 0.1400438 5e-07 -0.0001322 0 -0.0000107 -6.55e-05 0.0000205 NA -2.9e-06 NA NA NA 9 14318.98 -28619.94 23.30585 8.693194e-06
1 0.1403646 6e-07 -0.0001456 0 -0.0000101 -6.41e-05 -0.0000496 0e+00 NA NA NA NA 9 14314.57 -28611.12 32.12707 1.056029e-07
4 0.1406812 5e-07 -0.0001498 0 -0.0000090 -6.46e-05 -0.0000932 NA NA 0 NA NA 9 14311.61 -28605.18 38.06112 5.433915e-09
24 0.1403752 4e-07 -0.0001429 0 -0.0000912 -3.20e-05 -0.0000638 NA NA NA 1.57e-05 -4.9e-06 10 14309.73 -28599.43 43.81620 3.057838e-10
10 0.1398115 4e-07 -0.0001218 0 -0.0000886 -5.83e-05 0.0000503 NA -4.2e-06 NA 1.54e-05 NA 10 14309.45 -28598.87 44.37770 2.309328e-10
18 0.1400307 5e-07 -0.0001317 0 -0.0000116 -5.49e-05 0.0000244 NA -3.1e-06 NA NA -2.2e-06 10 14308.60 -28597.17 46.07873 9.865323e-11
9 0.1401851 8e-07 -0.0001399 0 -0.0000977 -5.53e-05 -0.0000401 -1e-07 NA NA 1.74e-05 NA 10 14305.24 -28590.45 52.79353 3.435659e-12
17 0.1403727 6e-07 -0.0001459 0 -0.0000104 -6.01e-05 -0.0000505 0e+00 NA NA NA -8.0e-07 10 14304.22 -28588.41 54.83621 1.237228e-12
3 0.1400829 5e-07 -0.0001345 0 -0.0000106 -6.52e-05 0.0000102 0e+00 -2.4e-06 NA NA NA 10 14303.99 -28587.95 55.29497 9.836251e-13
12 0.1406457 4e-07 -0.0001460 0 -0.0000867 -5.68e-05 -0.0001020 NA NA 0 1.54e-05 NA 10 14302.08 -28584.13 59.11472 1.456738e-13
20 0.1406925 5e-07 -0.0001502 0 -0.0000097 -5.45e-05 -0.0000943 NA NA 0 NA -2.1e-06 10 14301.22 -28582.40 60.84777 6.124294e-14
6 0.1404446 5e-07 -0.0001413 0 -0.0000094 -6.54e-05 -0.0000452 NA -1.7e-06 0 NA NA 10 14300.94 -28581.84 61.40234 4.641208e-14
26 0.1397606 4e-07 -0.0001200 0 -0.0000963 -2.97e-05 0.0000627 NA -4.8e-06 NA 1.65e-05 -5.8e-06 11 14299.18 -28576.32 66.92478 2.933925e-15
5 0.1406271 6e-07 -0.0001493 0 -0.0000090 -6.46e-05 -0.0000857 0e+00 NA 0 NA NA 10 14296.54 -28573.04 70.19997 5.704935e-16
25 0.1402101 8e-07 -0.0001405 0 -0.0001000 -4.18e-05 -0.0000431 -1e-07 NA NA 1.76e-05 -2.7e-06 11 14294.92 -28567.79 75.45124 4.130008e-17
11 0.1399502 7e-07 -0.0001306 0 -0.0000978 -5.63e-05 0.0000099 -1e-07 -2.0e-06 NA 1.73e-05 NA 11 14294.66 -28567.26 75.97923 3.171760e-17
19 0.1400554 5e-07 -0.0001332 0 -0.0000114 -5.66e-05 0.0000178 0e+00 -2.8e-06 NA NA -1.8e-06 11 14293.69 -28565.33 77.91047 1.207637e-17
28 0.1406735 4e-07 -0.0001469 0 -0.0000936 -2.97e-05 -0.0001054 NA NA 0 1.64e-05 -5.5e-06 11 14291.80 -28561.55 81.69801 1.817536e-18
14 0.1402437 4e-07 -0.0001315 0 -0.0000891 -5.80e-05 -0.0000209 NA -3.0e-06 0 1.58e-05 NA 11 14291.45 -28560.86 82.38049 1.292064e-18
22 0.1404325 5e-07 -0.0001409 0 -0.0000103 -5.41e-05 -0.0000414 NA -1.9e-06 0 NA -2.3e-06 11 14290.56 -28559.07 84.17102 5.278068e-19
13 0.1404731 8e-07 -0.0001438 0 -0.0000991 -5.56e-05 -0.0000802 -1e-07 NA 0 1.79e-05 NA 11 14287.28 -28552.51 90.73045 1.986613e-20
21 0.1406418 6e-07 -0.0001496 0 -0.0000094 -5.83e-05 -0.0000874 0e+00 NA 0 NA -1.3e-06 11 14286.19 -28550.34 92.90097 6.711053e-21
7 0.1404987 6e-07 -0.0001444 0 -0.0000092 -6.51e-05 -0.0000590 0e+00 -1.0e-06 0 NA NA 11 14285.96 -28549.87 93.37864 5.285250e-21
27 0.1398904 7e-07 -0.0001277 0 -0.0001010 -3.82e-05 0.0000257 -1e-07 -2.8e-06 NA 1.76e-05 -3.7e-06 12 14284.39 -28544.72 98.52532 4.031609e-22
30 0.1401981 4e-07 -0.0001298 0 -0.0000971 -2.82e-05 -0.0000092 NA -3.5e-06 0 1.69e-05 -6.0e-06 12 14281.20 -28538.34 104.90764 1.657971e-23
29 0.1405088 8e-07 -0.0001446 0 -0.0001020 -3.91e-05 -0.0000847 -1e-07 NA 0 1.82e-05 -3.4e-06 12 14276.97 -28529.88 113.36307 2.418259e-25
15 0.1404189 8e-07 -0.0001417 0 -0.0000991 -5.59e-05 -0.0000689 -1e-07 -4.0e-07 0 1.79e-05 NA 12 14276.69 -28529.33 113.91713 1.833119e-25
23 0.1404714 5e-07 -0.0001431 0 -0.0000099 -5.66e-05 -0.0000515 0e+00 -1.4e-06 0 NA -1.8e-06 12 14275.65 -28527.25 115.99526 6.485305e-26
31 0.1403591 7e-07 -0.0001388 0 -0.0001024 -3.76e-05 -0.0000531 -1e-07 -1.2e-06 0 1.82e-05 -3.8e-06 13 14266.42 -28506.78 136.46229 2.331148e-30

From the dredged model list we can see that the top model is one with no interactions. We ran this model on the MCC tree and 100 trees, noting that each model uses a unique set of values for \(\lambda_{DR}\)/\(\lambda_{ND}\) and a unique tree in the correlation structure.

#In both cases the top model is 1/2/3/4/5/6 no interaction terms. With no models within delta < 4: 

#Run model for DR
MCC.DR.top <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(0.985, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

saveRDS(MCC.DR.top, 'data/MCC.DR.top.rds')

#Run model for ND
MCC.ND.top <- gls(MCC.ND ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corPagel(0.9996, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

saveRDS(MCC.ND.top, 'data/MCC.ND.top.rds')
#Run the 100 models for DR and ND using the best model:


# No longer used now that we use 1000 trees

# #Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
# DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
#   left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
#             x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ES")), 
#             by = "TipLabel")
# })
# 
# #PGLS needs tiplabel as rowname
# DR.model.data <- lapply(DR.model.data, function(x) {
#   tibble::column_to_rownames(x, "TipLabel")})
# 
# #Prune the trees
# pruned.trees<-lapply(passerine.trees, function(x) {
#   drop.tip(x,x$tip.label[-match(restricted.data$TipLabel, x$tip.label)])
# })
# 
# #Use mapply to create a list of PGLS global models
# DR.pgls.models <- mcmapply(function(x,y) {
#   gls(ES ~ SDi 
#          + log(range.size.m2)
#          + bioclim4 #Seasonal variation
#          + residuals.PC1 #Spatial variation
#          + PC1.LIG #Long-term climate variation
#          + NPP,
#     corPagel(0.985, phy = y, fixed = TRUE), 
#     data = x, 
#     method = "REML")
# }, x = DR.model.data, y = pruned.trees,
# SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
# mc.cores = 8) #Specify core number 
# 
# saveRDS(DR.pgls.models, "data/DR.pgls.models.rds")
# 
# #Now for Node Density:
# ND.model.data <- lapply(nd.list, function(x) {
#   left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
#             x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
#             by = "TipLabel")
# })
# 
# #PGLS needs tiplabel as rowname
# ND.model.data <- lapply(ND.model.data, function(x) {
#   tibble::column_to_rownames(x, "TipLabel")})
# 
# #Use mapply to create a list of PGLS global models
# ND.pgls.models <- mcmapply(function(x,y) {
# gls(ND ~ SDi 
#          + log(range.size.m2)
#          + bioclim4 #Seasonal variation
#          + residuals.PC1 #Spatial variation
#          + PC1.LIG #Long-term climate variation
#          + NPP,
#     corPagel(0.9996, phy = y, fixed = TRUE), 
#     data = x, 
#     method = "REML")
# }, x = ND.model.data, y = pruned.trees, 
# SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
# mc.cores = 8) #Specify core number 
# 
# saveRDS(ND.pgls.models, "data/ND.pgls.models.rds")

From these models we obtained a distribution of the estimates for the models alongside the 95 % CIs of the MCC model. This enabled us to show the variation between trees, and the variation associated with the top MCC model. Given the variability across trees we ran 1000 ND and DR models, each allowed to estimate a phylogenetic signal (\(\lambda\)). The models were saved as an .rds object and read back into R with key coefficients extracted. Below we provide the code for this:

##########################################################
####### Read in DR models and extract estimates: #########
##########################################################

#We ran DR and ND models on 1000 trees
files.DR.SD <- list.files(path = "/Users/justincally/Dropbox/Runs Spartan/DR_SD/", pattern = "\\.rds$", full.names = TRUE) #1000 models
df <- list()
lapply(files.DR.SD, function(x) {
  tryCatch({
  model.x <- readRDS(x)
  name.x <- str_sub(x, -19, -5)
  df[[name.x]] <<- data.frame(model.x$coefficients,
                              confint(model.x),
                              coef(summary(model.x))[,2], #Std.Error
                              coef(summary(model.x))[,3], #t-val
                              coef(summary(model.x))[,4], #pval
                              model.x[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = name.x) %>% tibble::rownames_to_column()
  rm(model.x)
  gc(verbose = FALSE)
  },
  error = function(e) NULL
  )
})
DR.pgls.summary <- bind_rows(df) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name"))
#saveRDS(DR.pgls.summary, "data/DR.pgls.summary.rds") #Save a simple df with the model and coeff

##########################################################
####### Read in ND models and extract estimates: #########
##########################################################

files.ND.SD <- list.files(path = "/Users/justincally/Dropbox/Runs Spartan/ND_SD/", pattern = "\\.rds$", full.names = TRUE) #1000 models
df <- list()
lapply(files.ND.SD, function(x) {
  tryCatch({
  model.x <- readRDS(x) 
  name.x <- str_sub(x, -19, -5)
  df[[name.x]] <<- data.frame(model.x$coefficients, 
                              confint(model.x), 
                              coef(summary(model.x))[,2], #Std.Error
                              coef(summary(model.x))[,3], #t-val
                              coef(summary(model.x))[,4], #pval
                              model.x[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = name.x) %>% tibble::rownames_to_column()
  rm(model.x)
  gc(verbose = FALSE)
  },
  error = function(e) NULL
  )
})
ND.pgls.summary <- bind_rows(df) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name"))
#saveRDS(ND.pgls.summary, "data/ND.pgls.summary.rds") #Save a simple df with the model and coeff
#Read in the 1000 tree summary df and the MCC tree model
DR.pgls.summary <- readRDS("data/DR.pgls.summary.rds")
MCC.DR.top <- readRDS('data/MCC.DR.top.rds')
MCC.DR.summary <-  data.frame(MCC.DR.top$coefficients,
                              confint(MCC.DR.top),
                              coef(summary(MCC.DR.top))[,2], #Std.Error
                              coef(summary(MCC.DR.top))[,3], #t-val
                              coef(summary(MCC.DR.top))[,4], #pval
                              MCC.DR.top[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = "MCC_model") %>% 
  tibble::rownames_to_column() %>% 
  `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name"))

#Now plot this for Node Density (ND):  

ND.pgls.summary <- readRDS("data/ND.pgls.summary.rds")
MCC.ND.top <- readRDS('data/MCC.ND.top.rds')
MCC.ND.summary <- data.frame(MCC.ND.top$coefficients,
                              confint(MCC.ND.top),
                              coef(summary(MCC.ND.top))[,2], #Std.Error
                              coef(summary(MCC.ND.top))[,3], #t-val
                              coef(summary(MCC.ND.top))[,4], #pval
                              MCC.ND.top[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = "MCC_model") %>% 
  tibble::rownames_to_column() %>% 
  `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name"))

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `SDi` = "Sexual Dichromatism"
                    )



DR.plot <-DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  #geom_errorbarh(aes(xmin = LCI, xmax = UCI, colour = Parameter), position = position_jitter(seed = 1), height = 0)+
  geom_point(shape = 21, alpha = 0.5, size = 0.75, position = position_jitter(seed = 1))+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

DR.plot <- DR.plot + geom_errorbarh(data = MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("DR Models")


ND.plot <-ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

ND.plot <- ND.plot + geom_errorbarh(data = MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("ND Models")

PGLS Models on BAMM estimates

Using the previously generated BAMM dataframe we undertook model selection for speciation and extinction. The process is similar to that used on DR and ND, except given the use of a Bayesian approach in BAMM we can make use of varying levels of uncertainty between tips (species) by constructing a weighted model, where the weight is the inverse of the variance, such that more precise estimates of speciation or extinction at a given tip (species) holds higher weight in the model.

Based on preliminary findings we found that Pagel’s lambda was = 1 and running corPagel lead to problems of convergence. Therefore we ran the following models assuming Brownian Motion with corBrownian.

MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')

#Create model dataframe for use in models
MCC.BAMM.model.data <- left_join(restricted.data %>% 
                                   dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
                                 MCC.BAMM.df %>% as.data.frame(), 
                                 by = "TipLabel")

saveRDS(MCC.BAMM.model.data, 'data/MCC.BAMM.model.data.rds')
#Prune tree
pruned.MCC.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(MCC.BAMM.model.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(MCC.BAMM.model.data) <- MCC.BAMM.model.data$TipLabel

#Run model for DR
MCC.BAMM.lambda <- gls(mean.lambda ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                weights = ~ sqrt(var.lambda), #sqrt to account for overdispersedskewed variance distribution
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")

#Run model for DR
MCC.BAMM.mu <- gls(mean.mu ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")


#Dredge the global MCC models

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("MCC.BAMM.model.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

dredged.MCC.lambda <- pdredge(MCC.BAMM.lambda, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(dredged.MCC.lambda, "data/dredged.MCC.lambda.rds")

dredged.MCC.mu <- pdredge(MCC.BAMM.mu, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(dredged.MCC.mu, "data/dredged.MCC.mu.rds")

Table S6: The dredged models for BAMM speciation and BAMM extinction both show the top model is one with no interactions, with \(\delta AICc > 4\). This is the same situation as DR/ND models (see above).

dredged.MCC.lambda <- readRDS("data/dredged.MCC.lambda.rds")
dredged.MCC.mu <- readRDS("data/dredged.MCC.mu.rds")
kable(dredged.MCC.lambda, "html", caption = "BAMM-Speciation Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM-Speciation Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 SDi bioclim4:SDi log(range.size.m2):SDi NPP:SDi PC1.LIG:SDi residuals.PC1:SDi df logLik AICc delta weight
0 -2.109224 -9.0e-07 -0.0002226 0 -0.0000486 0.0000032 -0.0000143 NA NA NA NA NA 8 6820.146 -13624.27 0.00000 9.998913e-01
16 -2.109173 -9.0e-07 -0.0002238 0 -0.0000506 0.0000665 -0.0000198 NA NA NA NA -1.36e-05 9 6811.008 -13603.98 20.28182 3.942868e-05
8 -2.109262 -1.0e-06 -0.0002208 0 -0.0000967 0.0000102 -0.0000161 NA NA NA 1.04e-05 NA 9 6810.969 -13603.91 20.36033 3.791075e-05
2 -2.109394 -9.0e-07 -0.0002161 0 -0.0000491 0.0000029 0.0000283 NA -1.6e-06 NA NA NA 9 6810.768 -13603.50 20.76241 3.100651e-05
1 -2.109086 -1.2e-06 -0.0002252 0 -0.0000509 0.0000034 -0.0000299 1e-07 NA NA NA NA 9 6806.362 -13594.69 29.57314 3.786400e-07
4 -2.108965 -9.0e-07 -0.0002265 0 -0.0000477 0.0000071 -0.0000513 NA NA 0 NA NA 9 6803.212 -13588.39 35.87361 1.622171e-08
24 -2.109214 -1.0e-06 -0.0002217 0 -0.0001182 0.0000930 -0.0000238 NA NA NA 1.45e-05 -1.72e-05 10 6801.890 -13583.74 40.52471 1.585342e-09
18 -2.109479 -9.0e-07 -0.0002121 0 -0.0000516 0.0000677 0.0000568 NA -2.9e-06 NA NA -1.40e-05 10 6801.639 -13583.24 41.02738 1.233023e-09
10 -2.109492 -1.0e-06 -0.0002120 0 -0.0000981 0.0000099 0.0000414 NA -2.2e-06 NA 1.06e-05 NA 10 6801.594 -13583.15 41.11749 1.178701e-09
17 -2.108914 -1.3e-06 -0.0002287 0 -0.0000553 0.0000805 -0.0000494 1e-07 NA NA NA -1.65e-05 10 6797.286 -13574.53 49.73208 1.587699e-11
9 -2.109189 -1.1e-06 -0.0002223 0 -0.0000941 0.0000098 -0.0000239 0e+00 NA NA 9.60e-06 NA 10 6797.221 -13574.40 49.86347 1.486751e-11
3 -2.109504 -1.3e-06 -0.0002081 0 -0.0000533 0.0000026 0.0000821 1e-07 -4.5e-06 NA NA NA 10 6797.092 -13574.15 50.12130 1.306924e-11
20 -2.108935 -9.0e-07 -0.0002274 0 -0.0000498 0.0000685 -0.0000540 NA NA 0 NA -1.32e-05 10 6794.072 -13568.11 56.15975 6.382903e-13
12 -2.108990 -1.0e-06 -0.0002249 0 -0.0000975 0.0000145 -0.0000553 NA NA 0 1.08e-05 NA 10 6794.038 -13568.04 56.22814 6.168331e-13
6 -2.108997 -9.0e-07 -0.0002254 0 -0.0000478 0.0000070 -0.0000436 NA -3.0e-07 0 NA NA 10 6793.853 -13567.67 56.59952 5.122981e-13
26 -2.109643 -1.0e-06 -0.0002052 0 -0.0001217 0.0000956 0.0000833 NA -4.0e-06 NA 1.50e-05 -1.78e-05 11 6792.528 -13563.01 61.25528 4.995015e-14
5 -2.108810 -1.2e-06 -0.0002294 0 -0.0000502 0.0000074 -0.0000692 1e-07 NA 0 NA NA 10 6789.431 -13558.82 65.44261 6.155593e-15
25 -2.109027 -1.3e-06 -0.0002253 0 -0.0001138 0.0000998 -0.0000441 1e-07 NA NA 1.29e-05 -1.88e-05 11 6788.180 -13554.32 69.95146 6.459320e-16
19 -2.109728 -1.5e-06 -0.0001949 0 -0.0000608 0.0000916 0.0001737 2e-07 -9.0e-06 NA NA -1.92e-05 11 6788.073 -13554.10 70.16681 5.799929e-16
11 -2.109549 -1.2e-06 -0.0002075 0 -0.0000948 0.0000089 0.0000734 1e-07 -3.9e-06 NA 9.30e-06 NA 11 6787.951 -13553.86 70.41085 5.133695e-16
28 -2.108962 -1.0e-06 -0.0002254 0 -0.0001186 0.0000956 -0.0000602 NA NA 0 1.48e-05 -1.69e-05 11 6784.957 -13547.87 76.39784 2.572603e-17
22 -2.109125 -9.0e-07 -0.0002205 0 -0.0000504 0.0000691 -0.0000078 NA -1.7e-06 0 NA -1.35e-05 11 6784.722 -13547.40 76.86884 2.032805e-17
14 -2.109081 -1.0e-06 -0.0002215 0 -0.0000980 0.0000143 -0.0000332 NA -8.0e-07 0 1.08e-05 NA 11 6784.681 -13547.32 76.95100 1.950988e-17
21 -2.108658 -1.4e-06 -0.0002326 0 -0.0000546 0.0000829 -0.0000859 1e-07 NA 0 NA -1.62e-05 11 6780.353 -13538.66 85.60628 2.575051e-19
13 -2.108906 -1.2e-06 -0.0002265 0 -0.0000946 0.0000141 -0.0000644 0e+00 NA 0 9.90e-06 NA 11 6780.291 -13538.54 85.72956 2.421109e-19
7 -2.109117 -1.3e-06 -0.0002174 0 -0.0000519 0.0000066 0.0000104 1e-07 -3.0e-06 0 NA NA 11 6780.176 -13538.31 85.95995 2.157677e-19
27 -2.109816 -1.5e-06 -0.0001925 0 -0.0001182 0.0001103 0.0001722 1e-07 -8.7e-06 NA 1.27e-05 -2.14e-05 12 6778.966 -13533.88 90.38940 2.355832e-20
30 -2.109280 -1.0e-06 -0.0002137 0 -0.0001210 0.0000971 0.0000171 NA -2.8e-06 0 1.51e-05 -1.74e-05 12 6775.613 -13527.17 97.09545 8.240321e-22
29 -2.108763 -1.3e-06 -0.0002293 0 -0.0001140 0.0001027 -0.0000819 1e-07 NA 0 1.31e-05 -1.86e-05 12 6771.249 -13518.44 105.82259 1.049228e-23
23 -2.109414 -1.5e-06 -0.0002026 0 -0.0000595 0.0000920 0.0001139 2e-07 -7.7e-06 0 NA -1.86e-05 12 6771.155 -13518.26 106.01095 9.549222e-24
15 -2.109143 -1.2e-06 -0.0002173 0 -0.0000950 0.0000133 -0.0000026 0e+00 -2.4e-06 0 9.70e-06 NA 12 6771.038 -13518.02 106.24408 8.498514e-24
31 -2.109484 -1.4e-06 -0.0002008 0 -0.0001177 0.0001110 0.0001088 1e-07 -7.4e-06 0 1.29e-05 -2.08e-05 13 6762.050 -13498.04 126.23036 3.884888e-28
kable(dredged.MCC.mu, "html", caption = "BAMM-Extinction Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM-Extinction Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 SDi bioclim4:SDi log(range.size.m2):SDi NPP:SDi PC1.LIG:SDi residuals.PC1:SDi df logLik AICc delta weight
0 -4.355734 -1.3e-06 -0.0004223 0e+00 -0.0002270 -0.0000772 0.0000238 NA NA NA NA NA 8 2587.287 -5158.550 0.00000 9.997938e-01
16 -4.355867 -1.2e-06 -0.0004235 0e+00 -0.0002362 0.0000601 0.0000108 NA NA NA NA -2.95e-05 9 2578.874 -5139.717 18.83277 8.136278e-05
8 -4.355674 -1.2e-06 -0.0004250 0e+00 -0.0001830 -0.0000845 0.0000266 NA NA NA -8.60e-06 NA 9 2578.673 -5139.316 19.23424 6.656529e-05
2 -4.356956 -1.2e-06 -0.0003781 0e+00 -0.0002330 -0.0000786 0.0002791 NA -9.60e-06 NA NA NA 9 2578.528 -5139.024 19.52580 5.753574e-05
1 -4.355547 -1.7e-06 -0.0004254 0e+00 -0.0002278 -0.0000776 0.0000003 1e-07 NA NA NA NA 9 2574.172 -5130.313 28.23706 7.384306e-07
4 -4.354968 -1.2e-06 -0.0004346 -1e-07 -0.0002228 -0.0000758 -0.0000587 NA NA 0 NA NA 9 2570.981 -5123.932 34.61805 3.038754e-08
24 -4.355832 -1.2e-06 -0.0004249 0e+00 -0.0002130 0.0000523 0.0000126 NA NA NA -4.50e-06 -2.87e-05 10 2570.269 -5120.500 38.05030 5.462512e-09
18 -4.357425 -1.1e-06 -0.0003676 0e+00 -0.0002444 0.0000667 0.0003337 NA -1.21e-05 NA NA -3.13e-05 10 2570.128 -5120.217 38.33279 4.742985e-09
10 -4.356860 -1.2e-06 -0.0003821 0e+00 -0.0001912 -0.0000855 0.0002737 NA -9.30e-06 NA -8.10e-06 NA 10 2569.914 -5119.790 38.76023 3.830315e-09
17 -4.355488 -2.0e-06 -0.0004305 0e+00 -0.0002398 0.0000863 -0.0000426 2e-07 NA NA NA -3.53e-05 10 2565.824 -5111.611 46.93907 6.415370e-11
9 -4.355411 -1.8e-06 -0.0004298 0e+00 -0.0001722 -0.0000871 -0.0000035 1e-07 NA NA -1.09e-05 NA 10 2565.588 -5111.138 47.41236 5.063458e-11
3 -4.357274 -2.0e-06 -0.0003596 0e+00 -0.0002379 -0.0000803 0.0003756 2e-07 -1.49e-05 NA NA NA 10 2565.518 -5110.997 47.55267 4.720417e-11
20 -4.355073 -1.1e-06 -0.0004363 -1e-07 -0.0002321 0.0000643 -0.0000753 NA NA 0 NA -3.01e-05 10 2562.571 -5105.104 53.44577 2.479179e-12
12 -4.354931 -1.2e-06 -0.0004367 -1e-07 -0.0001835 -0.0000824 -0.0000544 NA NA 0 -7.70e-06 NA 10 2562.367 -5104.697 53.85297 2.022492e-12
6 -4.355908 -1.2e-06 -0.0004020 -1e-07 -0.0002275 -0.0000769 0.0001309 NA -6.80e-06 0 NA NA 10 2562.235 -5104.433 54.11716 1.772220e-12
26 -4.357372 -1.1e-06 -0.0003696 0e+00 -0.0002256 0.0000603 0.0003301 NA -1.20e-05 NA -3.60e-06 -3.06e-05 11 2561.524 -5101.002 57.54826 3.187608e-13
5 -4.354731 -1.7e-06 -0.0004385 -1e-07 -0.0002236 -0.0000762 -0.0000880 1e-07 NA 0 NA NA 10 2557.870 -5095.701 62.84860 2.251693e-14
19 -4.358259 -2.6e-06 -0.0003256 0e+00 -0.0002588 0.0001202 0.0005522 3e-07 -2.40e-05 NA NA -4.36e-05 11 2557.241 -5092.437 66.11295 4.402161e-15
25 -4.355393 -2.1e-06 -0.0004335 0e+00 -0.0001996 0.0000750 -0.0000442 2e-07 NA NA -7.80e-06 -3.44e-05 11 2557.241 -5092.436 66.11435 4.399069e-15
11 -4.357198 -2.1e-06 -0.0003617 0e+00 -0.0001778 -0.0000907 0.0003868 2e-07 -1.55e-05 NA -1.18e-05 NA 11 2556.937 -5091.828 66.72150 3.247290e-15
28 -4.355054 -1.1e-06 -0.0004372 -1e-07 -0.0002144 0.0000583 -0.0000730 NA NA 0 -3.40e-06 -2.94e-05 11 2553.967 -5085.888 72.66211 1.665462e-16
22 -4.356368 -1.1e-06 -0.0003916 -1e-07 -0.0002388 0.0000689 0.0001843 NA -9.40e-06 0 NA -3.14e-05 11 2553.836 -5085.626 72.92376 1.461222e-16
14 -4.355842 -1.2e-06 -0.0004051 -1e-07 -0.0001893 -0.0000833 0.0001290 NA -6.60e-06 0 -7.40e-06 NA 11 2553.622 -5085.198 73.35223 1.179443e-16
21 -4.354607 -2.1e-06 -0.0004447 -1e-07 -0.0002357 0.0000929 -0.0001392 2e-07 NA 0 NA -3.64e-05 11 2549.528 -5077.011 81.53875 1.967864e-18
13 -4.354622 -1.8e-06 -0.0004423 -1e-07 -0.0001718 -0.0000851 -0.0000897 1e-07 NA 0 -1.01e-05 NA 11 2549.285 -5076.524 82.02569 1.542617e-18
7 -4.356238 -2.0e-06 -0.0003834 -1e-07 -0.0002323 -0.0000786 0.0002277 2e-07 -1.21e-05 0 NA NA 11 2549.225 -5076.404 82.14595 1.452594e-18
27 -4.358181 -2.7e-06 -0.0003279 0e+00 -0.0002149 0.0001082 0.0005563 4e-07 -2.42e-05 NA -8.50e-06 -4.26e-05 12 2548.659 -5073.264 85.28577 3.022330e-19
30 -4.356335 -1.1e-06 -0.0003930 -1e-07 -0.0002240 0.0000638 0.0001826 NA -9.30e-06 0 -2.80e-06 -3.08e-05 12 2545.233 -5066.411 92.13852 9.823906e-21
23 -4.357230 -2.6e-06 -0.0003493 -1e-07 -0.0002533 0.0001216 0.0004053 3e-07 -2.12e-05 0 NA -4.35e-05 12 2540.948 -5057.842 100.70759 1.353739e-22
29 -4.354536 -2.1e-06 -0.0004472 -1e-07 -0.0002001 0.0000828 -0.0001391 2e-07 NA 0 -6.90e-06 -3.56e-05 12 2540.944 -5057.835 100.71495 1.348764e-22
15 -4.356204 -2.0e-06 -0.0003845 -1e-07 -0.0001765 -0.0000883 0.0002434 2e-07 -1.27e-05 0 -1.10e-05 NA 12 2540.644 -5057.234 101.31603 9.986502e-23
31 -4.357186 -2.6e-06 -0.0003508 -1e-07 -0.0002136 0.0001106 0.0004127 4e-07 -2.15e-05 0 -7.70e-06 -4.26e-05 13 2532.366 -5038.669 119.88067 9.292931e-27

Given that we ran models weighted according to the variance of the response (\(\lambda_{BAMM}\) and \(\mu_{BAMM}\)), we checked to see whether a weighted model is favourable to an unweighted model using an anova to compare AIC values.

#Run the top model for the MCC
#Run model for ND
MCC.Lambda.top <- gls(mean.lambda ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Lambda.top, 'data/MCC.Lambda.top.rds')

MCC.Mu.top <- gls(mean.mu ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")

saveRDS(MCC.Mu.top, 'data/MCC.Mu.top.rds')


#We can also see how the models look without the weightings: 
MCC.Lambda.top.unweighted <- gls(mean.lambda ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Lambda.top.unweighted, 'data/MCC.Lambda.top.unweighted.rds')


MCC.Mu.top.unweighted <- gls(mean.mu ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Mu.top.unweighted, 'data/MCC.Mu.top.unweighted.rds')

Table S6: Models with a set of weights for the response based on the inverse of the variance of the posterior distribution in \(\lambda_{BAMM}\) and \(\mu_{BAMM}\) have much lower AIC values, indicating that the weighting scheme improves the fit of the model. In any case, we found that an unweighted model did not change qualitative results.

MCC.Lambda.top <- readRDS('data/MCC.Lambda.top.rds')
MCC.Lambda.top.unweighted <- readRDS('data/MCC.Lambda.top.unweighted.rds')
MCC.Mu.top <- readRDS('data/MCC.Mu.top.rds')
MCC.Mu.top.unweighted <- readRDS('data/MCC.Mu.top.unweighted.rds')

anova(MCC.Lambda.top, MCC.Lambda.top.unweighted) %>% pander(split.table = Inf)
  call Model df AIC BIC logLik
MCC.Lambda.top gls(model = mean.lambda ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), weights = ~sqrt(var.lambda), method = “REML”) 1 8 -13624 -13571 6820
MCC.Lambda.top.unweighted gls(model = mean.lambda ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), method = “REML”) 2 8 -11788 -11735 5902
anova(MCC.Mu.top, MCC.Mu.top.unweighted) %>% pander(split.table = Inf)
  call Model df AIC BIC logLik
MCC.Mu.top gls(model = mean.mu ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), weights = ~sqrt(var.mu), method = “REML”) 1 8 -5159 -5105 2587
MCC.Mu.top.unweighted gls(model = mean.mu ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), method = “REML”) 2 8 -5320 -5266 2668

We ran the top (no interactions) model on the 100 trees, each with unique estimates of speciation and extinction from the BAMM runs.

#Read in the BAMM data for the 100 trees
BAMM.df <- readRDS('data/BAMM.df.rds')

#Take the restricted data and make it simpler with just responses and predictors.Note that we join the BAMM for the 100 trees

BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
BAMM.model.data <- lapply(BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(restricted.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(mean.lambda ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.lambda),      
    corBrownian(phy = y), 
    data = x, 
    method = "REML")
}, x = BAMM.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(BAMM.lambda.pgls.models, "data/BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(mean.mu ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.mu),      
    corBrownian(phy = y), 
    data = x, 
    method = "REML")
}, x = BAMM.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(BAMM.mu.pgls.models, "data/BAMM.mu.pgls.models.rds")
BAMM.lambda.pgls.models <- readRDS("data/BAMM.lambda.pgls.models.rds")

BAMM.lambda.pgls.summary <- lapply(BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

#Now plot this for Lambda

MCC.Lambda.top <- readRDS('data/MCC.Lambda.top.rds')
MCC.lambda.summary <- data.frame(MCC.Lambda.top$coefficients,
                              confint(MCC.Lambda.top),
                              coef(summary(MCC.Lambda.top))[,2], #Std.Error
                              coef(summary(MCC.Lambda.top))[,3], #t-val
                              coef(summary(MCC.Lambda.top))[,4], #pval
                              MCC.Lambda.top[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = "MCC_model") %>% 
  tibble::rownames_to_column() %>% 
  `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name"))


parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `SDi` = "Sexual Dichromatism"
                    )

BAMM.lambda.pgls.summary <- bind_rows(BAMM.lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

  # for (x in unique(BAMM.lambda.pgls.summary$Parameter)[2:7]){
  #   filter(Parameter == x & between(Estimate, left = as.numeric(hpd.Lambda.top[1,x]), right = as.numeric(hpd.Lambda.top[2,x])))
  #   } 


remove_outliers <- function(x, na.rm = TRUE, ...) {
  qnt <- quantile(x, probs=c(.25, .75), na.rm = na.rm, ...)
  H <- 1.5 * IQR(x, na.rm = na.rm)
  y <- x
  y[x < (qnt[1] - H)] <- NA
  y[x > (qnt[2] + H)] <- NA
  y
}

BAMM.lambda.pgls.summary.RO <- dcast(BAMM.lambda.pgls.summary %>% filter(Parameter != "(Intercept)"), Estimate ~ Parameter, value.var = "Estimate")
BAMM.lambda.pgls.summary.RO$Estimate <- NULL
BAMM.lambda.pgls.summary.RO <- sapply(BAMM.lambda.pgls.summary.RO, function(x) {
  remove_outliers(x, na.rm = T)})
BAMM.lambda.pgls.summary.RO <-melt(BAMM.lambda.pgls.summary.RO) %>% na.omit()
BAMM.lambda.pgls.summary.RO$Var1 <- NULL
colnames(BAMM.lambda.pgls.summary.RO) <- c("Parameter", "Estimate")

BAMM.lambda.plot <- BAMM.lambda.pgls.summary.RO %>%
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

BAMM.lambda.plot <- BAMM.lambda.plot + geom_errorbarh(data = MCC.lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 7, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("BAMM Speciation")


#Figure for extinction

BAMM.mu.pgls.models <- readRDS("data/BAMM.mu.pgls.models.rds")

BAMM.mu.pgls.summary <- lapply(BAMM.mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

BAMM.mu.pgls.summary <- bind_rows(BAMM.mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

#Now plot this for Lambda

MCC.Mu.top <- readRDS('data/MCC.Mu.top.rds')
MCC.mu.summary <- data.frame(MCC.Mu.top$coefficients,
                              confint(MCC.Mu.top),
                              coef(summary(MCC.Mu.top))[,2], #Std.Error
                              coef(summary(MCC.Mu.top))[,3], #t-val
                              coef(summary(MCC.Mu.top))[,4], #pval
                              MCC.Mu.top[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = "MCC_model") %>% 
  tibble::rownames_to_column() %>% 
  `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name"))

#Get rid of outliers
BAMM.mu.pgls.summary.RO <- dcast(BAMM.mu.pgls.summary %>% select(Parameter,Estimate) %>% filter(Parameter != "(Intercept)"), Estimate ~ Parameter, value.var = "Estimate")
BAMM.mu.pgls.summary.RO$Estimate <- NULL
BAMM.mu.pgls.summary.RO <- sapply(BAMM.mu.pgls.summary.RO, function(x) {
  remove_outliers(x, na.rm = T)})
BAMM.mu.pgls.summary.RO <-melt(BAMM.mu.pgls.summary.RO) %>% na.omit()
BAMM.mu.pgls.summary.RO$Var1 <- NULL
colnames(BAMM.mu.pgls.summary.RO) <- c("Parameter", "Estimate")

BAMM.mu.plot <-BAMM.mu.pgls.summary.RO %>%
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

BAMM.mu.plot <- BAMM.mu.plot + geom_errorbarh(data = MCC.mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 7, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,0.5,.5),"cm"))+
  ggtitle("BAMM Extinction")

PGLS Model Summary

grid.arrange(
  symmetrise_scale(DR.plot, "x"),
  symmetrise_scale(ND.plot, "x"),
  symmetrise_scale(BAMM.lambda.plot, "x"),
  symmetrise_scale(BAMM.mu.plot, "x"),
  nrow = 1
) 

Figure S9: This figure is the same basic figure as seen in the manuscript (Figure 1). It provides model estimates for four response variables across 100 random trees alongside the MCC tree.

Table S7: The estimates of the MCC models plotted above are based on the following data tables.

MCC.DR.summary %>% select(-model_name) %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC DR Estimates", split.table = Inf)
MCC DR Estimates
Parameter Estimate LCI UCI SE tval pval lambda
SDi -0.001279 -0.003324 0.0007664 0.001043 -1.226 0.2204 0.985
log(range.size.m2) -0.006579 -0.01078 -0.002381 0.002142 -3.072 0.002136 0.985
bioclim4 2.137e-05 -3.369e-05 7.644e-05 2.809e-05 0.7607 0.4469 0.985
residuals.PC1 -0.002441 -0.007567 0.002686 0.002616 -0.9332 0.3507 0.985
PC1.LIG 0.001592 -0.003672 0.006855 0.002685 0.5927 0.5534 0.985
NPP 2.902e-07 -1.325e-06 1.905e-06 8.239e-07 0.3522 0.7247 0.985
MCC.ND.summary %>% select(-model_name) %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC ND Estimates", split.table = Inf)
MCC ND Estimates
Parameter Estimate LCI UCI SE tval pval lambda
SDi -5.745e-05 -0.0001239 8.955e-06 3.388e-05 -1.696 0.09 0.9996
log(range.size.m2) -0.0001462 -0.0002826 -9.713e-06 6.963e-05 -2.099 0.03582 0.9996
bioclim4 4.488e-07 -1.428e-06 2.325e-06 9.575e-07 0.4688 0.6393 0.9996
residuals.PC1 -6.402e-05 -0.0002359 0.0001079 8.769e-05 -0.7301 0.4654 0.9996
PC1.LIG -1.016e-05 -0.0001763 0.000156 8.479e-05 -0.1198 0.9046 0.9996
NPP -9.174e-09 -6.464e-08 4.629e-08 2.83e-08 -0.3242 0.7458 0.9996
MCC.lambda.summary %>% select(-model_name) %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Speciation Estimates", split.table = Inf)
MCC BAMM Speciation Estimates
Parameter Estimate LCI UCI SE tval pval lambda
SDi -1.429e-05 -0.0002574 0.0002288 0.000124 -0.1152 0.9083 1
log(range.size.m2) -0.0002226 -0.0006901 0.0002448 0.0002385 -0.9335 0.3506 1
bioclim4 -9.384e-07 -7.233e-06 5.356e-06 3.212e-06 -0.2922 0.7702 1
residuals.PC1 3.219e-06 -0.0005798 0.0005862 0.0002975 0.01082 0.9914 1
PC1.LIG -4.856e-05 -0.0005749 0.0004778 0.0002685 -0.1808 0.8565 1
NPP -2.088e-08 -2.282e-07 1.864e-07 1.058e-07 -0.1974 0.8435 1
MCC.mu.summary %>% select(-model_name) %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Extinction Estimates", split.table = Inf)
MCC BAMM Extinction Estimates
Parameter Estimate LCI UCI SE tval pval lambda
SDi 2.385e-05 -0.0004099 0.0004575 0.0002213 0.1078 0.9142 1
log(range.size.m2) -0.0004223 -0.00133 0.0004853 0.0004631 -0.9119 0.3619 1
bioclim4 -1.283e-06 -1.399e-05 1.142e-05 6.482e-06 -0.1979 0.8431 1
residuals.PC1 -7.718e-05 -0.001233 0.001078 0.0005896 -0.1309 0.8958 1
PC1.LIG -0.000227 -0.001284 0.0008298 0.0005392 -0.421 0.6737 1
NPP -2.036e-08 -4.179e-07 3.772e-07 2.028e-07 -0.1004 0.9201 1

From the estimates of 100 trees we calculated Highest posterior density (HPD) intervals. Note that these do no account for the CI of each estimate. Thus it is a more a representation of tree uncertainty than model uncertainty.

Table S8: HPD intervals calculated for the above figure (i.e. models using RGB values of sexual dichromatism). The HPD range is determined using the hdi function with a 95 % credible interval. These intervals do not take into account the variance associated with each interval and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that is in the same direction (+ or -).

hpd.DR.top <- list()
for(x in unique(DR.pgls.summary$Parameter)) {
hpd.DR.top[[x]] = hdi(DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select("Estimate"))
}
hpd.DR.top <- bind_rows(hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)")


saveRDS(hpd.DR.top, 'data/hpd.DR.top.rds')

#For ND
hpd.ND.top <- list()
for(x in unique(ND.pgls.summary$Parameter)) {
hpd.ND.top[[x]] = hdi(ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.ND.top <- bind_rows(hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)") 

saveRDS(hpd.ND.top, 'data/hpd.ND.top.rds')

hpd.Lambda.top <- list()
for(x in unique(BAMM.lambda.pgls.summary$Parameter)) {
hpd.Lambda.top[[x]] = hdi(BAMM.lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.Lambda.top <- bind_rows(hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)") 


saveRDS(hpd.Lambda.top, 'data/hpd.Lambda.top.rds')

hpd.Mu.top <- list()
for(x in unique(BAMM.mu.pgls.summary$Parameter)) {
hpd.Mu.top[[x]] = hdi(BAMM.mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.Mu.top <- bind_rows(hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)")


saveRDS(hpd.Mu.top, 'data/hpd.Mu.top.rds')

hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "DR HPD Interval")
DR HPD Interval
  SDi log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.00193 -0.00844 -2.14e-05 -0.00329 -0.00388 -1.4e-06
Upper 0.00157 -0.00182 5.62e-05 0.0048 0.00524 1.53e-06
hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "ND HPD Interval")
ND HPD Interval
  SDi log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -7.68e-05 -0.00019 -1.13e-06 -9.11e-05 -0.00015 -4.26e-08
Upper 5.97e-05 9.12e-06 1.46e-06 0.000175 0.000117 4.28e-08
hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "BAMM Speciation HPD Interval")
BAMM Speciation HPD Interval
  SDi log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.00831 -0.00785 -0.000178 -0.0146 -0.0106 -5.54e-06
Upper 0.00474 0.0115 0.000116 0.0185 0.00835 2.43e-06
hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "BAMM Extinction HPD Interval")
BAMM Extinction HPD Interval
  SDi log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.00788 -0.0166 -0.000158 -0.0216 -0.0161 -3.96e-06
Upper 0.0121 0.024 0.000423 0.0368 0.0263 9.84e-06

From the HPD 95 % intervals and the intervals for the 95 % CI we can gain an indication of how much variability there is within a model in comparison to variability between trees. To do this we calculated the width of the HPD interval relative to the 95 % CI interval from the MCC tree.

Table S9: Model estimates from \(\lambda_{BAMM}\) and \(\mu_{BAMM}\) tip-rate estimates show high levels of variability between trees relative to estimates from \(\lambda_{DR}\) and \(\lambda_{ND}\). The interval ratio is calculated as the width of the HPD 95 % CI relative to the MCC 95 % CI. A value of 1 suggests that model estimates from 95 % of the trees fall within the MCC 95 % CI. Here, the large values for the estimates from BAMM tip-rates suggest that there is high variation in tip-rate estimates across trees that is not see within the DR and ND models.

hpd.DR.top2 <- hpd.DR.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.DR.top2$HPD.Interval <- hpd.DR.top2$Upper - hpd.DR.top2$Lower
DR.intervals <- left_join(hpd.DR.top2, MCC.DR.summary, by = "Parameter")
DR.intervals$MCC.Interval <- DR.intervals$UCI - DR.intervals$LCI
DR.intervals$DR.IntervalRatio <- hpd.DR.top2$HPD.Interval/DR.intervals$MCC.Interval

hpd.ND.top2 <- hpd.ND.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.ND.top2$HPD.Interval <- hpd.ND.top2$Upper - hpd.ND.top2$Lower
ND.intervals <- left_join(hpd.ND.top2, MCC.ND.summary, by = "Parameter")
ND.intervals$MCC.Interval <- ND.intervals$UCI - ND.intervals$LCI
ND.intervals$ND.IntervalRatio <- hpd.ND.top2$HPD.Interval/ND.intervals$MCC.Interval


hpd.Lambda.top2 <- hpd.Lambda.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.Lambda.top2$HPD.Interval <- hpd.Lambda.top2$Upper - hpd.Lambda.top2$Lower
Lambda.intervals <- left_join(hpd.Lambda.top2, MCC.lambda.summary, by = "Parameter")
Lambda.intervals$MCC.Interval <- Lambda.intervals$UCI - Lambda.intervals$LCI
Lambda.intervals$Lambda.IntervalRatio <- hpd.Lambda.top2$HPD.Interval/Lambda.intervals$MCC.Interval


hpd.Mu.top2 <- hpd.Mu.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.Mu.top2$HPD.Interval <- hpd.Mu.top2$Upper - hpd.Mu.top2$Lower
Mu.intervals <- left_join(hpd.Mu.top2, MCC.mu.summary, by = "Parameter")
Mu.intervals$MCC.Interval <- Mu.intervals$UCI - Mu.intervals$LCI
Mu.intervals$Mu.IntervalRatio <- hpd.Mu.top2$HPD.Interval/Mu.intervals$MCC.Interval


plyr::join_all(list(DR.intervals %>% select(Parameter, DR.IntervalRatio),
          ND.intervals %>% select(Parameter, ND.IntervalRatio),
          Lambda.intervals %>% select(Parameter, Lambda.IntervalRatio),
          Mu.intervals %>% select(Parameter, Mu.IntervalRatio)), by = "Parameter", type = "left") %>% `colnames<-`(c("Parameter", "λDR.Interval.Ratio", "λND.Interval.Ratio", "λBAMM.Interval.Ratio", "μBAMM.Interval.Ratio")) %>%
  pander(digits = 3, split.table = Inf) 
Parameter λDR.Interval.Ratio λND.Interval.Ratio λBAMM.Interval.Ratio μBAMM.Interval.Ratio
SDi 0.856 1.03 26.9 23
log(range.size.m2) 0.788 0.729 20.7 22.3
bioclim4 0.705 0.689 23.3 22.9
residuals.PC1 0.79 0.774 28.3 25.3
PC1.LIG 0.867 0.805 18 20.1
NPP 0.906 0.77 19.2 17.4

Subsetted analysis with spectrophotometry data

Using the dataset from Armenta et al. (2008) we conducted the analysis on a subset of species for which dichromatism values were derived from spectrophotometry data. The measure of dichromatism is a difference measure between the sexes, based on a model of bird colour vision. To make this dataset comparable with the RGB measure of dichromatism, we use the absolute difference between the sexes; thereby making the scale from monochromatism to dichromatism rather than female colouration to male colouration.

#Read in Armenta data
Armenta.data <- read.csv('data/Armenta_2008.csv', stringsAsFactors = F)
MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')
#A couple of rows have "-", remove from dataset
Armenta.data <- Armenta.data %>% dplyr::select(binomial, Colour.discriminability) %>% mutate(Colour.discriminability = replace(Colour.discriminability, Colour.discriminability == "-", "NA")) %>% filter(Colour.discriminability != "NA")
Armenta.data$Colour.discriminability <- Armenta.data$Colour.discriminability %>% as.numeric()

MCC.Armenta.model.data <- inner_join(restricted.data %>% dplyr::select(TipLabel, binomial, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP, MCC.DR, MCC.ND), Armenta.data, by = "binomial")

MCC.Armenta.model.data <- inner_join(MCC.Armenta.model.data, MCC.BAMM.df, by = "TipLabel") %>% filter(Colour.discriminability != "NA")

#Plot the correlation
p1 <- MCC.Armenta.model.data %>% ggplot(aes(x = SDi, y = Colour.discriminability))+
  geom_point()+
  geom_smooth(method = 'loess')+
  theme_minimal()

grid.newpage()
ggExtra::ggMarginal(
  p = p1,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = 'gray'
)

Figure S10: There is a correlation between the RGB measures and the Spectophotometry measures. The RGB measures seem to be more noisy around the lower values.

R-squared and correlation for relationship between RGB and Spectrophotometry Data:

data_frame(
R2 = summary(lm(SDi ~ Colour.discriminability,
   data = MCC.Armenta.model.data))$r.squared,
r = cor(MCC.Armenta.model.data$SDi, MCC.Armenta.model.data$Colour.discriminability)) %>% pander()
R2 r
0.6166 0.7853

Run PGLS model for the data

#Prune tree
pruned.MCC.Armenta.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(MCC.Armenta.model.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(MCC.Armenta.model.data) <- MCC.Armenta.model.data$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.Armenta.global <- gls(MCC.DR ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = FALSE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.ND.Armenta.global <- gls(MCC.ND ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), #lambda = 1
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Lambda.Armenta.global <- gls(mean.lambda ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Extinction.Armenta.global <- gls(mean.mu ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                weights = ~ sqrt(var.mu),
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), 
                data = MCC.Armenta.model.data, 
                method = "REML")
#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("MCC.Armenta.model.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.Armenta.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

Armenta.dredged.ND.model <- pdredge(MCC.ND.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.ND.model, 'data/Armenta.dredged.ND.model.rds')

Armenta.dredged.DR.model <- pdredge(MCC.DR.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.DR.model, 'data/Armenta.dredged.DR.model.rds')

Armenta.dredged.spec.model <- pdredge(MCC.Lambda.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.spec.model, 'data/Armenta.dredged.spec.model.rds')

Armenta.dredged.extinct.model <- pdredge(MCC.Extinction.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.extinct.model, 'data/Armenta.dredged.extinct.model.rds')

Table S10: The dredged models using a subset of species for which sexual dichromatism was measured using spectrophotometry. All four show the top model is one with no interactions, with \(\delta AICc > 4\).

Armenta.dredged.ND.model <- readRDS("data/Armenta.dredged.ND.model.rds")
Armenta.dredged.DR.model <- readRDS("data/Armenta.dredged.DR.model.rds")
Armenta.dredged.spec.model <- readRDS("data/Armenta.dredged.spec.model.rds")
Armenta.dredged.extinct.model <- readRDS("data/Armenta.dredged.extinct.model.rds")


kable(Armenta.dredged.ND.model, "html", caption = "ND Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
ND Spectrophotometry Dichromatism Dredge Table
(Intercept) bioclim4 Colour.discriminability log(range.size.m2) NPP PC1.LIG residuals.PC1 bioclim4:Colour.discriminability Colour.discriminability:log(range.size.m2) Colour.discriminability:NPP Colour.discriminability:PC1.LIG Colour.discriminability:residuals.PC1 df logLik AICc delta weight
0 0.2078007 7.0e-06 0.0000823 -0.0000392 5e-07 -0.0008872 0.0010554 NA NA NA NA NA 8 1033.2993 -2050.347 0.00000 9.972277e-01
8 0.2129124 6.3e-06 0.0003663 -0.0002059 4e-07 -0.0008461 0.0010488 NA NA NA -0.0009035 NA 9 1027.8922 -2037.469 12.87782 1.593716e-03
2 0.2060950 7.3e-06 -0.0081443 0.0000163 5e-07 -0.0008871 0.0010736 NA 0.0002873 NA NA NA 9 1027.0518 -2035.788 14.55847 6.878007e-04
16 0.2084122 6.9e-06 0.0001265 -0.0000615 5e-07 -0.0008772 0.0010554 NA NA NA NA -0.0001579 9 1026.6976 -2035.080 15.26687 4.826543e-04
1 0.2072802 7.1e-06 -0.0002068 -0.0000215 5e-07 -0.0008912 0.0010556 6.00e-07 NA NA NA NA 9 1021.9507 -2025.586 24.76079 4.188489e-06
10 0.2094495 7.1e-06 -0.0258624 -0.0000932 5e-07 -0.0008298 0.0011044 NA 0.0009198 NA -0.0012567 NA 10 1022.3449 -2024.304 26.04315 2.205954e-06
24 0.2127459 6.3e-06 0.0003407 -0.0001967 4e-07 -0.0008557 0.0010478 NA NA NA -0.0010482 0.0002544 10 1021.4636 -2022.541 27.80565 9.138503e-07
4 0.2077667 7.2e-06 0.0008705 -0.0000400 5e-07 -0.0009178 0.0010595 NA NA -1e-07 NA NA 9 1019.5545 -2020.794 29.55320 3.814161e-07
18 0.2067050 7.2e-06 -0.0081183 -0.0000059 5e-07 -0.0008770 0.0010736 NA 0.0002879 NA NA -0.0001586 10 1020.4512 -2020.516 29.83052 3.320313e-07
9 0.2096091 7.9e-06 -0.0054190 -0.0000819 4e-07 -0.0008690 0.0010431 1.38e-05 NA NA -0.0023003 NA 10 1018.8088 -2017.232 33.11531 6.425347e-08
3 0.2062391 7.3e-06 -0.0088321 0.0000111 5e-07 -0.0008846 0.0010754 -4.00e-07 0.0003176 NA NA NA 10 1015.8402 -2011.295 39.05241 3.301204e-09
17 0.2074343 7.2e-06 -0.0005586 -0.0000289 5e-07 -0.0008821 0.0010559 1.60e-06 NA NA NA -0.0002338 10 1015.4920 -2010.598 39.74893 2.330361e-09
26 0.2085390 7.3e-06 -0.0305168 -0.0000567 5e-07 -0.0008443 0.0011122 NA 0.0010814 NA -0.0015814 0.0004620 11 1016.1463 -2009.829 40.51838 1.586137e-09
12 0.2134049 6.8e-06 0.0028561 -0.0002277 4e-07 -0.0009367 0.0010606 NA NA -3e-07 -0.0010093 NA 10 1014.3515 -2008.317 42.02978 7.449790e-10
6 0.2060863 7.3e-06 -0.0082061 0.0000166 5e-07 -0.0008865 0.0010736 NA 0.0002889 0e+00 NA NA 10 1013.4408 -2006.496 43.85126 2.996502e-10
20 0.2083967 7.1e-06 0.0009743 -0.0000630 5e-07 -0.0009097 0.0010598 NA NA -1e-07 NA -0.0001634 10 1012.9578 -2005.530 44.81722 1.848667e-10
11 0.2079243 8.2e-06 -0.0200882 -0.0000287 4e-07 -0.0008572 0.0010759 1.24e-05 0.0005345 NA -0.0023670 NA 11 1012.8361 -2003.208 47.13875 5.790884e-11
25 0.2095828 7.9e-06 -0.0053570 -0.0000797 4e-07 -0.0008726 0.0010427 1.36e-05 NA NA -0.0023417 0.0001034 11 1012.3085 -2002.153 48.19395 3.416727e-11
19 0.2066033 7.3e-06 -0.0071140 -0.0000023 5e-07 -0.0008789 0.0010710 6.00e-07 0.0002440 NA NA -0.0001866 11 1009.3880 -1996.312 54.03494 1.841860e-12
5 0.2075081 7.2e-06 0.0006229 -0.0000309 5e-07 -0.0009159 0.0010590 3.00e-07 NA -1e-07 NA NA 10 1008.2776 -1996.169 54.17762 1.715038e-12
14 0.2096741 7.1e-06 -0.0243702 -0.0001015 5e-07 -0.0008457 0.0011041 NA 0.0008821 0e+00 -0.0012600 NA 11 1008.7374 -1995.011 55.33603 9.610122e-13
28 0.2132517 6.9e-06 0.0030119 -0.0002188 4e-07 -0.0009543 0.0010603 NA NA -3e-07 -0.0011822 0.0002903 11 1007.9550 -1993.446 56.90089 4.394636e-13
22 0.2067366 7.2e-06 -0.0079025 -0.0000071 5e-07 -0.0008793 0.0010736 NA 0.0002824 0e+00 NA -0.0001589 11 1006.8425 -1991.221 59.12593 1.444645e-13
27 0.2075208 8.3e-06 -0.0229037 -0.0000127 4e-07 -0.0008635 0.0010817 1.17e-05 0.0006425 NA -0.0024795 0.0002476 12 1006.4516 -1988.354 61.99296 3.445044e-14
13 0.2097077 8.0e-06 -0.0050441 -0.0000860 4e-07 -0.0008800 0.0010446 1.36e-05 NA 0e+00 -0.0022976 NA 11 1005.1214 -1987.779 62.56802 2.584171e-14
7 0.2062628 7.3e-06 -0.0086970 0.0000103 5e-07 -0.0008861 0.0010754 -4.00e-07 0.0003144 0e+00 NA NA 11 1002.2442 -1982.024 68.32243 1.454681e-15
21 0.2075790 7.2e-06 0.0000033 -0.0000348 5e-07 -0.0008988 0.0010582 1.30e-06 NA 0e+00 NA -0.0002248 11 1001.8222 -1981.180 69.16645 9.538725e-16
30 0.2087086 7.3e-06 -0.0293939 -0.0000630 5e-07 -0.0008561 0.0011120 NA 0.0010529 0e+00 -0.0015830 0.0004607 12 1002.5378 -1980.526 69.82065 6.877521e-16
15 0.2073599 8.1e-06 -0.0234613 -0.0000077 4e-07 -0.0008205 0.0010759 1.27e-05 0.0006139 1e-07 -0.0023855 NA 12 999.2623 -1973.975 76.37154 2.599714e-17
29 0.2097143 8.0e-06 -0.0048443 -0.0000851 4e-07 -0.0008877 0.0010448 1.34e-05 NA 0e+00 -0.0023413 0.0001115 12 998.6337 -1972.718 77.62878 1.386499e-17
23 0.2065974 7.3e-06 -0.0071475 -0.0000021 5e-07 -0.0008785 0.0010710 6.00e-07 0.0002448 0e+00 NA -0.0001866 12 995.7933 -1967.037 83.30953 8.097675e-19
31 0.2069695 8.1e-06 -0.0261965 0.0000078 4e-07 -0.0008274 0.0010817 1.20e-05 0.0007198 1e-07 -0.0024970 0.0002463 13 992.8778 -1959.114 91.23332 1.540747e-20
kable(Armenta.dredged.DR.model, "html", caption = "DR Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
DR Spectrophotometry Dichromatism Dredge Table
(Intercept) bioclim4 Colour.discriminability log(range.size.m2) NPP PC1.LIG residuals.PC1 bioclim4:Colour.discriminability Colour.discriminability:log(range.size.m2) Colour.discriminability:NPP Colour.discriminability:PC1.LIG Colour.discriminability:residuals.PC1 df logLik AICc delta weight
0 -2.677449 0.0001793 0.0070637 0.0090357 7.8e-06 -0.0316263 0.0157111 NA NA NA NA NA 9 -512.3713 1043.058 0.000000 9.479264e-01
8 -2.631082 0.0001609 0.0092736 0.0078425 7.2e-06 -0.0289217 0.0168116 NA NA NA -0.0169400 NA 10 -514.7191 1049.824 6.766463 3.217035e-02
2 -2.669646 0.0001774 0.0667327 0.0087989 7.8e-06 -0.0315029 0.0156588 NA -0.0020914 NA NA NA 10 -515.8283 1052.042 8.984711 1.061132e-02
16 -2.676194 0.0001781 0.0078263 0.0089900 7.9e-06 -0.0313609 0.0158720 NA NA NA NA -0.0018026 10 -516.0825 1052.551 9.493220 8.229020e-03
10 -2.653267 0.0001645 -0.2225684 0.0085929 7.4e-06 -0.0288817 0.0172761 NA 0.0081403 NA -0.0199814 NA 11 -517.9109 1058.286 15.228039 4.677867e-04
24 -2.624764 0.0001607 0.0065382 0.0077914 6.9e-06 -0.0293425 0.0164764 NA NA NA -0.0213142 0.0079015 11 -518.0505 1058.565 15.507264 4.068319e-04
18 -2.668515 0.0001762 0.0666983 0.0087573 7.8e-06 -0.0312406 0.0158197 NA -0.0020637 NA NA -0.0017903 11 -519.5385 1061.541 18.483252 9.187289e-05
1 -2.672231 0.0001758 0.0136928 0.0089013 7.8e-06 -0.0313185 0.0158244 -0.0000153 NA NA NA NA 10 -520.8853 1062.157 19.098835 6.753289e-05
9 -2.642694 0.0001831 -0.0713938 0.0081969 7.2e-06 -0.0295885 0.0167470 0.0001928 NA NA -0.0366782 NA 11 -521.5607 1065.585 22.527606 1.216094e-05
4 -2.674122 0.0001870 0.0486739 0.0088590 7.7e-06 -0.0332248 0.0161648 NA NA -4.6e-06 NA NA 10 -523.1226 1066.631 23.573445 7.208851e-06
26 -2.654292 0.0001659 -0.3226718 0.0088332 7.1e-06 -0.0294234 0.0170182 NA 0.0115291 NA -0.0269724 0.0103438 12 -521.0552 1066.660 23.601946 7.106850e-06
3 -2.669429 0.0001755 0.0432560 0.0088168 7.7e-06 -0.0313221 0.0157801 -0.0000117 -0.0010917 NA NA NA 11 -524.1766 1070.817 27.759395 8.889935e-07
17 -2.672734 0.0001757 0.0125002 0.0089099 7.8e-06 -0.0312094 0.0159148 -0.0000114 NA NA NA -0.0012551 11 -524.5169 1071.498 28.440062 6.325482e-07
12 -2.618684 0.0001704 0.0766293 0.0073732 6.9e-06 -0.0310819 0.0177078 NA NA -7.4e-06 -0.0195332 NA 11 -525.1095 1072.683 29.625169 3.497445e-07
11 -2.643318 0.0001829 -0.0777098 0.0082272 7.2e-06 -0.0295756 0.0167773 0.0001918 0.0002365 NA -0.0366666 NA 12 -524.8610 1074.271 31.213435 1.580751e-07
25 -2.636787 0.0001826 -0.0727548 0.0081488 6.9e-06 -0.0299694 0.0164410 0.0001900 NA NA -0.0404521 0.0073388 12 -524.9295 1074.408 31.350431 1.476098e-07
6 -2.638645 0.0001833 0.3287355 0.0077133 7.3e-06 -0.0336174 0.0161532 NA -0.0089933 -7.2e-06 NA NA 11 -526.2428 1074.950 31.891755 1.126077e-07
20 -2.673128 0.0001858 0.0488452 0.0088252 7.7e-06 -0.0329910 0.0162928 NA NA -4.6e-06 NA -0.0014566 11 -526.8371 1076.138 33.080334 6.215432e-08
19 -2.668577 0.0001755 0.0575069 0.0087750 7.8e-06 -0.0311954 0.0158572 -0.0000046 -0.0016749 NA NA -0.0015781 12 -527.7668 1080.083 37.025073 8.647332e-09
28 -2.607912 0.0001719 0.0851395 0.0072194 6.4e-06 -0.0320431 0.0174150 NA NA -8.8e-06 -0.0259111 0.0106718 12 -528.2282 1081.006 37.947974 5.451007e-09
14 -2.622800 0.0001705 0.0363708 0.0075223 6.9e-06 -0.0309746 0.0177577 NA 0.0013123 -7.1e-06 -0.0199138 NA 12 -528.3370 1081.223 38.165451 4.889364e-09
27 -2.644978 0.0001827 -0.1663156 0.0084481 7.0e-06 -0.0299489 0.0166173 0.0001785 0.0034451 NA -0.0409777 0.0080944 13 -528.1585 1082.959 39.901175 2.052790e-09
22 -2.638168 0.0001824 0.3265876 0.0076986 7.3e-06 -0.0334243 0.0162596 NA -0.0089199 -7.1e-06 NA -0.0011666 12 -529.9585 1084.466 41.408565 9.660930e-10
5 -2.656479 0.0001790 0.0857025 0.0083627 7.4e-06 -0.0328940 0.0166878 -0.0000479 NA -6.4e-06 NA NA 11 -531.4042 1085.272 42.214658 6.456216e-10
13 -2.635643 0.0001846 -0.0303830 0.0079411 7.1e-06 -0.0304933 0.0171656 0.0001682 NA -3.4e-06 -0.0353442 NA 12 -532.3308 1089.211 46.153047 9.010892e-11
30 -2.621005 0.0001725 -0.0485580 0.0076844 6.5e-06 -0.0317599 0.0175235 NA 0.0043689 -7.8e-06 -0.0275285 0.0112744 13 -531.4040 1089.450 46.392211 7.995286e-11
7 -2.635611 0.0001789 0.2868764 0.0076675 7.2e-06 -0.0333063 0.0165021 -0.0000309 -0.0068797 -7.8e-06 NA NA 12 -534.5255 1093.600 50.542593 1.003668e-11
21 -2.655184 0.0001787 0.0903183 0.0083517 7.4e-06 -0.0330515 0.0166707 -0.0000538 NA -6.7e-06 NA 0.0013286 12 -534.9923 1094.534 51.476130 6.293250e-12
15 -2.627334 0.0001844 0.0525129 0.0076658 7.0e-06 -0.0306837 0.0170951 0.0001727 -0.0027887 -4.0e-06 -0.0349672 NA 13 -535.5286 1097.699 54.641420 1.292829e-12
29 -2.625333 0.0001847 -0.0143206 0.0077683 6.7e-06 -0.0313522 0.0169675 0.0001542 NA -4.8e-06 -0.0393824 0.0089632 13 -535.5787 1097.799 54.741494 1.229731e-12
23 -2.635476 0.0001787 0.2850029 0.0076828 7.2e-06 -0.0333509 0.0165134 -0.0000337 -0.0067505 -7.9e-06 NA 0.0005011 13 -538.1025 1102.847 59.789133 9.856650e-14
31 -2.625599 0.0001845 -0.0171436 0.0077861 6.7e-06 -0.0313401 0.0169870 0.0001536 0.0001046 -4.8e-06 -0.0393741 0.0089753 14 -538.7648 1106.272 63.213836 1.778544e-14
kable(Armenta.dredged.spec.model, "html", caption = "BAMM Speciation Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM Speciation Spectrophotometry Dichromatism Dredge Table
(Intercept) bioclim4 Colour.discriminability log(range.size.m2) NPP PC1.LIG residuals.PC1 bioclim4:Colour.discriminability Colour.discriminability:log(range.size.m2) Colour.discriminability:NPP Colour.discriminability:PC1.LIG Colour.discriminability:residuals.PC1 df logLik AICc delta weight
0 -2.481836 -2.57e-05 -0.0005284 0.0018331 -1.6e-06 -0.0042427 0.0024700 NA NA NA NA NA 8 181.4716 -346.6914 0.00000 9.919978e-01
2 -2.487950 -2.40e-05 -0.0426377 0.0020014 -1.6e-06 -0.0042761 0.0025397 NA 0.0014635 NA NA NA 9 176.7417 -335.1681 11.52334 3.120684e-03
8 -2.477096 -2.70e-05 -0.0001201 0.0017016 -1.7e-06 -0.0041154 0.0025774 NA NA NA -0.0011822 NA 9 176.5728 -334.8304 11.86107 2.635792e-03
16 -2.478548 -2.62e-05 -0.0000828 0.0017268 -1.6e-06 -0.0040970 0.0025549 NA NA NA NA -0.0009388 9 176.3917 -334.4682 12.22323 2.199227e-03
1 -2.474805 -2.77e-05 0.0040932 0.0016238 -1.7e-06 -0.0041103 0.0025673 -1.04e-05 NA NA NA NA 9 171.7110 -325.1068 21.58465 2.039218e-05
10 -2.483988 -2.51e-05 -0.0748327 0.0018904 -1.6e-06 -0.0040672 0.0027923 NA 0.0026087 NA -0.0021814 NA 10 172.1199 -323.8538 22.83759 1.089913e-05
18 -2.484660 -2.46e-05 -0.0421396 0.0018950 -1.6e-06 -0.0041305 0.0026243 NA 0.0014617 NA NA -0.0009373 10 171.6621 -322.9382 23.75321 6.895528e-06
24 -2.476484 -2.69e-05 0.0000327 0.0016753 -1.7e-06 -0.0040635 0.0025976 NA NA NA -0.0008290 -0.0005790 10 171.5863 -322.7867 23.90470 6.392512e-06
4 -2.481930 -2.60e-05 -0.0016972 0.0018368 -1.6e-06 -0.0042004 0.0024630 NA NA 1e-07 NA NA 9 169.2266 -320.1380 26.55344 1.700215e-06
3 -2.480943 -2.62e-05 -0.0779854 0.0017784 -1.6e-06 -0.0040584 0.0027984 -1.98e-05 0.0029986 NA NA NA 10 167.3535 -314.3211 32.37035 9.276363e-08
9 -2.474764 -2.77e-05 0.0034705 0.0016249 -1.7e-06 -0.0041004 0.0025779 -8.80e-06 NA NA -0.0002837 NA 10 167.1595 -313.9330 32.75845 7.640185e-08
17 -2.474401 -2.76e-05 0.0034228 0.0016072 -1.7e-06 -0.0040570 0.0025943 -8.30e-06 NA NA NA -0.0005120 10 166.6820 -312.9780 33.71338 4.739625e-08
26 -2.484021 -2.51e-05 -0.0750209 0.0018915 -1.6e-06 -0.0040685 0.0027923 NA 0.0026151 NA -0.0021936 0.0000159 11 167.1583 -311.8527 34.83873 2.700082e-08
6 -2.490882 -2.50e-05 -0.0657857 0.0020885 -1.6e-06 -0.0040492 0.0025269 NA 0.0020376 7e-07 NA NA 10 164.6815 -308.9771 37.71432 6.411356e-09
12 -2.476878 -2.68e-05 0.0007949 0.0016947 -1.7e-06 -0.0041442 0.0025862 NA NA -1e-07 -0.0012184 NA 10 164.3628 -308.3397 38.35172 4.661658e-09
20 -2.478642 -2.65e-05 -0.0012702 0.0017306 -1.6e-06 -0.0040539 0.0025477 NA NA 1e-07 NA -0.0009392 10 164.1475 -307.9091 38.78232 3.758692e-09
11 -2.481001 -2.61e-05 -0.0817896 0.0017853 -1.6e-06 -0.0040314 0.0028320 -1.59e-05 0.0030789 NA -0.0007324 NA 11 162.8230 -303.1820 43.50945 3.536345e-10
19 -2.481122 -2.62e-05 -0.0790795 0.0017843 -1.6e-06 -0.0040690 0.0027961 -2.04e-05 0.0030438 NA NA 0.0001094 11 162.3497 -302.2354 44.45598 2.203009e-10
25 -2.474401 -2.76e-05 0.0033410 0.0016076 -1.7e-06 -0.0040563 0.0025955 -8.10e-06 NA NA -0.0000417 -0.0005045 11 162.1738 -301.8835 44.80788 1.847576e-10
5 -2.473758 -2.73e-05 0.0073858 0.0015905 -1.7e-06 -0.0041943 0.0025954 -1.16e-05 NA -3e-07 NA NA 10 159.5559 -298.7259 47.96553 3.810027e-11
14 -2.486651 -2.59e-05 -0.0943958 0.0019693 -1.6e-06 -0.0038730 0.0027753 NA 0.0030858 6e-07 -0.0021311 NA 11 160.0492 -297.6343 49.05709 2.207492e-11
22 -2.487594 -2.55e-05 -0.0653495 0.0019822 -1.6e-06 -0.0039028 0.0026116 NA 0.0020373 7e-07 NA -0.0009389 11 159.6029 -296.7418 49.94958 1.412855e-11
28 -2.476430 -2.69e-05 0.0002873 0.0016737 -1.7e-06 -0.0040721 0.0025999 NA NA 0e+00 -0.0008429 -0.0005729 11 159.3900 -296.3159 50.37548 1.141866e-11
27 -2.481516 -2.61e-05 -0.0857305 0.0018034 -1.6e-06 -0.0040549 0.0028333 -1.66e-05 0.0032240 NA -0.0009007 0.0003067 12 157.8890 -291.2286 55.46279 8.972611e-13
7 -2.482718 -2.66e-05 -0.0884878 0.0018317 -1.6e-06 -0.0039504 0.0027818 -1.90e-05 0.0032336 4e-07 NA NA 11 155.2726 -288.0811 58.61028 1.859732e-13
13 -2.473747 -2.73e-05 0.0067594 0.0015923 -1.7e-06 -0.0041836 0.0026040 -1.01e-05 NA -3e-07 -0.0002490 NA 11 155.0065 -287.5489 59.14249 1.425219e-13
21 -2.473679 -2.73e-05 0.0059146 0.0015847 -1.7e-06 -0.0041247 0.0026118 -9.50e-06 NA -2e-07 NA -0.0004535 11 154.5467 -286.6295 60.06196 8.999580e-14
30 -2.486613 -2.59e-05 -0.0941911 0.0019680 -1.6e-06 -0.0038712 0.0027753 NA 0.0030787 6e-07 -0.0021165 -0.0000190 12 155.0893 -285.6294 61.06205 5.458259e-14
15 -2.483019 -2.66e-05 -0.0941518 0.0018466 -1.6e-06 -0.0039058 0.0028173 -1.45e-05 0.0033548 4e-07 -0.0008224 NA 12 150.7529 -276.9566 69.73483 7.141405e-16
23 -2.482781 -2.66e-05 -0.0888749 0.0018339 -1.6e-06 -0.0039574 0.0027809 -1.93e-05 0.0032527 4e-07 NA 0.0000549 12 150.2768 -276.0043 70.68711 4.436087e-16
29 -2.473679 -2.73e-05 0.0058298 0.0015851 -1.7e-06 -0.0041240 0.0026131 -9.20e-06 NA -2e-07 -0.0000438 -0.0004456 12 150.0394 -275.5294 71.16202 3.498445e-16
31 -2.483365 -2.65e-05 -0.0969398 0.0018591 -1.6e-06 -0.0039314 0.0028190 -1.52e-05 0.0034656 4e-07 -0.0009611 0.0002601 13 145.8223 -265.0027 81.68872 1.811467e-18
kable(Armenta.dredged.extinct.model, "html", caption = "BAMM Extinction Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM Extinction Spectrophotometry Dichromatism Dredge Table
(Intercept) bioclim4 Colour.discriminability log(range.size.m2) NPP PC1.LIG residuals.PC1 bioclim4:Colour.discriminability Colour.discriminability:log(range.size.m2) Colour.discriminability:NPP Colour.discriminability:PC1.LIG Colour.discriminability:residuals.PC1 df logLik AICc delta weight
0 -5.637995 -2.50e-06 -0.0007061 0.0098071 -7.0e-07 -0.0158256 0.0047712 NA NA NA NA NA 8 -295.8228 607.8974 0.000000 9.760508e-01
8 -5.602968 -9.30e-06 0.0013779 0.0086260 -1.0e-06 -0.0153267 0.0048510 NA NA NA -0.0069576 NA 9 -299.3634 617.0421 9.144629 1.008652e-02
2 -5.660188 3.20e-06 -0.1229593 0.0104573 -5.0e-07 -0.0158528 0.0050603 NA 0.0042624 NA NA NA 9 -299.6657 617.6466 9.749147 7.455419e-03
16 -5.621902 -4.40e-06 0.0006272 0.0091847 -8.0e-07 -0.0155294 0.0049035 NA NA NA NA -0.0040005 9 -299.8668 618.0488 10.151407 6.097087e-03
10 -5.633627 -4.00e-07 -0.2645285 0.0094454 -8.0e-07 -0.0151327 0.0055227 NA 0.0093079 NA -0.0104941 NA 10 -302.5922 625.5703 17.672853 1.418604e-04
1 -5.607684 -1.07e-05 0.0181462 0.0086986 -8.0e-07 -0.0155623 0.0048618 -4.15e-05 NA NA NA NA 9 -304.4746 627.2645 19.367046 6.080954e-05
24 -5.600876 -9.20e-06 0.0016675 0.0085204 -1.0e-06 -0.0152675 0.0048977 NA NA NA -0.0059688 -0.0017575 10 -303.5391 627.4642 19.566765 5.503050e-05
18 -5.643552 1.10e-06 -0.1164458 0.0098227 -6.0e-07 -0.0155629 0.0051769 NA 0.0040807 NA NA -0.0039001 10 -303.7238 627.8335 19.936076 4.575184e-05
4 -5.638119 1.80e-06 0.0142059 0.0097805 -7.0e-07 -0.0164478 0.0048358 NA NA -1.6e-06 NA NA 9 -307.1995 632.7142 24.816808 3.986309e-06
3 -5.634607 -4.20e-06 -0.2350576 0.0093313 -6.0e-07 -0.0154353 0.0055561 -7.09e-05 0.0092936 NA NA NA 10 -307.7146 635.8152 27.917724 8.456988e-07
26 -5.634214 -3.00e-07 -0.2672034 0.0094678 -8.0e-07 -0.0151388 0.0055230 NA 0.0094002 NA -0.0106630 0.0002378 11 -306.7694 636.0028 28.105346 7.699705e-07
9 -5.601586 -9.90e-06 0.0039065 0.0085670 -1.0e-06 -0.0153312 0.0048572 -6.00e-06 NA NA -0.0063681 NA 10 -307.9862 636.3584 28.460957 6.445461e-07
17 -5.604788 -9.90e-06 0.0144636 0.0085713 -8.0e-07 -0.0154383 0.0049237 -3.16e-05 NA NA NA -0.0025285 10 -308.6159 637.6177 29.720240 3.434031e-07
12 -5.597067 -2.60e-06 0.0297843 0.0083691 -1.2e-06 -0.0164095 0.0049864 NA NA -3.1e-06 -0.0081761 NA 10 -310.5342 641.4544 33.557003 5.042685e-08
6 -5.658062 3.90e-06 -0.1065779 0.0103859 -5.0e-07 -0.0160366 0.0050515 NA 0.0038471 -5.0e-07 NA NA 10 -310.9655 642.3170 34.419559 3.276124e-08
20 -5.621692 2.00e-07 0.0167676 0.0091426 -8.0e-07 -0.0161953 0.0049761 NA NA -1.8e-06 NA -0.0040860 10 -311.2301 642.8462 34.948757 2.514469e-08
11 -5.629550 -2.60e-06 -0.2761375 0.0092269 -7.0e-07 -0.0151374 0.0056147 -2.96e-05 0.0101541 NA -0.0078861 NA 11 -311.1044 644.6727 36.775240 1.008859e-08
19 -5.633135 -4.20e-06 -0.2288248 0.0092821 -6.0e-07 -0.0154083 0.0055518 -6.77e-05 0.0090315 NA NA -0.0006231 11 -311.9005 646.2650 38.367539 4.550584e-09
25 -5.600250 -9.50e-06 0.0028834 0.0084940 -1.0e-06 -0.0152708 0.0048998 -2.90e-06 NA NA -0.0057034 -0.0017229 11 -312.1530 646.7699 38.872451 3.535305e-09
14 -5.629258 1.00e-06 -0.2337213 0.0092978 -9.0e-07 -0.0154856 0.0055093 NA 0.0085388 -1.0e-06 -0.0105809 NA 11 -313.8771 650.2182 42.320798 6.304168e-10
5 -5.597405 -4.60e-06 0.0563410 0.0082567 -1.0e-06 -0.0167909 0.0050304 -5.60e-05 NA -3.5e-06 NA NA 10 -315.5637 651.5134 43.616012 3.298950e-10
28 -5.595617 -2.70e-06 0.0290960 0.0082946 -1.1e-06 -0.0163281 0.0050187 NA NA -3.0e-06 -0.0073614 -0.0013782 11 -314.7225 651.9089 44.011480 2.707080e-10
22 -5.640154 2.10e-06 -0.0918395 0.0097077 -6.0e-07 -0.0158342 0.0051652 NA 0.0034565 -7.0e-07 NA -0.0039510 11 -315.0150 652.4939 44.596495 2.020536e-10
27 -5.631352 -2.50e-06 -0.2869084 0.0092903 -7.0e-07 -0.0151601 0.0056237 -3.21e-05 0.0105634 NA -0.0082829 0.0008701 12 -315.2429 655.0350 47.137607 5.671149e-11
7 -5.626340 -2.10e-06 -0.1850305 0.0090420 -7.0e-07 -0.0160338 0.0055485 -7.41e-05 0.0081220 -1.6e-06 NA NA 11 -318.9505 660.3650 52.467557 3.947181e-12
13 -5.591571 -3.90e-06 0.0418151 0.0081339 -1.2e-06 -0.0165464 0.0050234 -2.09e-05 NA -3.4e-06 -0.0062486 NA 11 -319.0833 660.6307 52.733234 3.456178e-12
30 -5.629761 1.00e-06 -0.2360819 0.0093169 -9.0e-07 -0.0154883 0.0055096 NA 0.0086185 -1.0e-06 -0.0107175 0.0001930 12 -318.0531 660.6554 52.758007 3.413633e-12
21 -5.595822 -4.40e-06 0.0510955 0.0081864 -1.0e-06 -0.0166179 0.0050674 -4.74e-05 NA -3.3e-06 NA -0.0019479 11 -319.7348 661.9335 54.036071 1.801725e-12
15 -5.622905 -9.00e-07 -0.2341699 0.0089926 -8.0e-07 -0.0156347 0.0056071 -3.32e-05 0.0091738 -1.3e-06 -0.0076925 NA 12 -322.3553 669.2600 61.362554 4.621288e-14
23 -5.625049 -2.10e-06 -0.1796751 0.0089989 -7.0e-07 -0.0160046 0.0055446 -7.10e-05 0.0078899 -1.6e-06 NA -0.0005723 12 -323.1361 670.8216 62.924154 2.116732e-14
29 -5.591062 -3.80e-06 0.0397984 0.0081020 -1.1e-06 -0.0164638 0.0050451 -1.83e-05 NA -3.3e-06 -0.0058219 -0.0011173 12 -323.2644 671.0781 63.180633 1.861968e-14
31 -5.624708 -9.00e-07 -0.2449379 0.0090561 -8.0e-07 -0.0156588 0.0056161 -3.57e-05 0.0095844 -1.3e-06 -0.0080923 0.0008777 13 -326.4929 679.6278 71.730326 2.590778e-16

Run the top model and 100 models for each one:

#In both cases the top model is 1/2/3/4/5/6 no interaction terms. With no models within delta < 4: 

#Run model for DR
Armenta.MCC.DR.top <- gls(MCC.DR ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = FALSE), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.DR.top, 'data/Armenta.MCC.DR.top.rds')

#Run model for ND
Armenta.MCC.ND.top <- gls(MCC.ND ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), #lambda = 1
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.ND.top, 'data/Armenta.MCC.ND.top.rds')

#Run the 100 models for DR and ND using the best model:
Armenta.data.noMCC <- inner_join(restricted.data %>% dplyr::select(TipLabel, binomial, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), Armenta.data, by = "binomial") %>% filter(Colour.discriminability != "NA")
#Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
Armenta.DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "DR")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.DR.model.data <- lapply(Armenta.DR.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
Armenta.pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(MCC.Armenta.model.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
Armenta.DR.pgls.models <- mcmapply(function(x,y) {
  gls(DR ~ Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    correlation = corPagel(1, phy = y, fixed = FALSE), 
    data = x, 
    method = "REML")
}, x = Armenta.DR.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.DR.pgls.models, "data/Armenta.DR.pgls.models.rds")

#Now for Node Density:
Armenta.ND.model.data <- lapply(nd.list, function(x) {
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.ND.model.data <- lapply(Armenta.ND.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Armenta.ND.pgls.models <- mcmapply(function(x,y) {
gls(ND ~ Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(1, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = Armenta.ND.model.data, y = Armenta.pruned.trees, 
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.ND.pgls.models, "data/Armenta.ND.pgls.models.rds")

#Run the BAMM models

Armenta.MCC.Lambda.top <- gls(mean.lambda ~ Colour.discriminability
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corBrownian(phy = pruned.MCC.Armenta.tree), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.Lambda.top, 'data/Armenta.MCC.Lambda.top.rds')

Armenta.MCC.Mu.top <- gls(mean.mu ~ Colour.discriminability
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.Armenta.tree), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.Mu.top, 'data/Armenta.MCC.Mu.top.rds')

Armenta.BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.BAMM.model.data <- lapply(Armenta.BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Armenta.BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(mean.lambda ~ Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.lambda),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Armenta.BAMM.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.BAMM.lambda.pgls.models, "data/Armenta.BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
Armenta.BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(mean.mu ~ Colour.discriminability
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.mu),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Armenta.BAMM.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.BAMM.mu.pgls.models, "data/Armenta.BAMM.mu.pgls.models.rds")

Using the MCC models as well as the 100 PGLS models we can generate the plots similar to Figure S8:

Armenta.DR.pgls.models <- readRDS("data/Armenta.DR.pgls.models.rds")
Armenta.MCC.DR.top <- readRDS('data/Armenta.MCC.DR.top.rds')

Armenta.DR.pgls.summary <- lapply(Armenta.DR.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.DR.summary <- data.frame(Armenta.MCC.DR.top$coefficients, confint(Armenta.MCC.DR.top)) %>% tibble::rownames_to_column()

Armenta.DR.pgls.summary <- bind_rows(Armenta.DR.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `Colour.discriminability` = "Sexual Dichromatism"
                    )

Armenta.DR.pgls.summary$Parameter = factor(Armenta.DR.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.DR.summary$Parameter = factor(Armenta.MCC.DR.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.DR.plot <-Armenta.DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.DR.plot <- Armenta.DR.plot + geom_errorbarh(data = Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta DR")

####################################
#For ND

Armenta.ND.pgls.models <- readRDS("data/Armenta.ND.pgls.models.rds")
Armenta.MCC.ND.top <- readRDS('data/Armenta.MCC.ND.top.rds')

Armenta.ND.pgls.summary <- lapply(Armenta.ND.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.ND.summary <- data.frame(Armenta.MCC.ND.top$coefficients, confint(Armenta.MCC.ND.top)) %>% tibble::rownames_to_column()

Armenta.ND.pgls.summary <- bind_rows(Armenta.ND.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.ND.pgls.summary$Parameter = factor(Armenta.ND.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.ND.summary$Parameter = factor(Armenta.MCC.ND.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.ND.plot <-Armenta.ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.ND.plot <- Armenta.ND.plot + geom_errorbarh(data = Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta ND")

######################

#For Lambda
Armenta.BAMM.lambda.pgls.models <- readRDS("data/Armenta.BAMM.lambda.pgls.models.rds")
Armenta.MCC.Lambda.top <- readRDS('data/Armenta.MCC.Lambda.top.rds')

Armenta.Lambda.pgls.summary <- lapply(Armenta.BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.Lambda.summary <- data.frame(Armenta.MCC.Lambda.top$coefficients, confint(Armenta.MCC.Lambda.top)) %>% tibble::rownames_to_column()

Armenta.Lambda.pgls.summary <- bind_rows(Armenta.Lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.Lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.Lambda.pgls.summary$Parameter = factor(Armenta.Lambda.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.Lambda.summary$Parameter = factor(Armenta.MCC.Lambda.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.Lambda.plot <-Armenta.Lambda.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.Lambda.plot <- Armenta.Lambda.plot + geom_errorbarh(data = Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta BAMM Speciaton")

######################

#For Mu

Armenta.BAMM.Mu.pgls.models <- readRDS("data/Armenta.BAMM.mu.pgls.models.rds")
Armenta.MCC.Mu.top <- readRDS('data/Armenta.MCC.Mu.top.rds')

Armenta.Mu.pgls.summary <- lapply(Armenta.BAMM.Mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.Mu.summary <- data.frame(Armenta.MCC.Mu.top$coefficients, confint(Armenta.MCC.Mu.top)) %>% tibble::rownames_to_column()

Armenta.Mu.pgls.summary <- bind_rows(Armenta.Mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.Mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.Mu.pgls.summary$Parameter = factor(Armenta.Mu.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.Mu.summary$Parameter = factor(Armenta.MCC.Mu.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.Mu.plot <-Armenta.Mu.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.Mu.plot <- Armenta.Mu.plot + geom_errorbarh(data = Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta BAMM Extinction")

Plot the results

grid.arrange(symmetrise_scale(Armenta.DR.plot, "x"),
             symmetrise_scale(Armenta.ND.plot, "x"),
             symmetrise_scale(Armenta.Lambda.plot, "x"), 
             symmetrise_scale(Armenta.Mu.plot, "x"), 
             nrow = 1)

Figure S11: Measures of sexual dichromatism using spectrophotometry do not change the main patterns. This dataset is a subset of the complete dataset (n = 581), thus drawing conclusions for the other predictors (e.g. borderline long term temperature variation and spatial temperature variation) potentially risks type I error. Additionally, these environmental predictors were not measured differently in this analysis compared to the analysis using RGB measures of dichromatism for all species so we refrain from drawing conclusions on the effects of the environmental predictors based on this smaller dataset.

Table S11: MCC estimates from the above plot are presented as numerical values below.

Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC DR Estimates from Spec measures")
MCC DR Estimates from Spec measures
Parameter Estimate LCI UCI
Colour.discriminability 0.007064 -0.04821 0.06233
log(range.size.m2) 0.009036 -0.02175 0.03982
bioclim4 0.0001793 -8.546e-05 0.0004441
residuals.PC1 0.01571 -0.007217 0.03864
PC1.LIG -0.03163 -0.07058 0.007327
NPP 7.831e-06 -5.329e-06 2.099e-05
Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC ND Estimates from Spec measures")
MCC ND Estimates from Spec measures
Parameter Estimate LCI UCI
Colour.discriminability 8.234e-05 -0.003517 0.003682
log(range.size.m2) -3.922e-05 -0.00174 0.001661
bioclim4 6.983e-06 -8.266e-06 2.223e-05
residuals.PC1 0.001055 -0.0001934 0.002304
PC1.LIG -0.0008872 -0.00306 0.001286
NPP 4.64e-07 -2.859e-07 1.214e-06
Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Speciation Estimates from Spec measures")
MCC BAMM Speciation Estimates from Spec measures
Parameter Estimate LCI UCI
Colour.discriminability -0.0005284 -0.01583 0.01478
log(range.size.m2) 0.001833 -0.005722 0.009388
bioclim4 -2.57e-05 -9.319e-05 4.179e-05
residuals.PC1 0.00247 -0.003289 0.008229
PC1.LIG -0.004243 -0.01356 0.005075
NPP -1.633e-06 -4.922e-06 1.657e-06
Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Extinction Estimates from Spec measures")
MCC BAMM Extinction Estimates from Spec measures
Parameter Estimate LCI UCI
Colour.discriminability -0.0007061 -0.03742 0.036
log(range.size.m2) 0.009807 -0.007571 0.02718
bioclim4 -2.453e-06 -0.0001574 0.0001525
residuals.PC1 0.004771 -0.008028 0.01757
PC1.LIG -0.01583 -0.03779 0.006138
NPP -6.892e-07 -8.33e-06 6.951e-06

Table S12: HPD intervals calculated for the above figure (i.e. models using Spectrophotometry values of sexual dichromatism). These intervals do not take into account the variance associated with each estimate and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that are in the same direction (+ or -). These intervals are calculated in the same way as in Table S7.

Armenta.hpd.DR.top <- list()
Armenta.DR.pgls.summary <- na.omit(Armenta.DR.pgls.summary)
for(x in unique(Armenta.DR.pgls.summary$Parameter)) {
Armenta.hpd.DR.top[[x]] = hdi(Armenta.DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.DR.top <- bind_rows(Armenta.hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper"))


saveRDS(Armenta.hpd.DR.top, 'data/Armenta.hpd.DR.top.rds')

#For ND
Armenta.hpd.ND.top <- list()
Armenta.ND.pgls.summary <- na.omit(Armenta.ND.pgls.summary)
for(x in unique(Armenta.ND.pgls.summary$Parameter)) {
Armenta.hpd.ND.top[[x]] = hdi(Armenta.ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.ND.top <- bind_rows(Armenta.hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Armenta.hpd.ND.top, 'data/Armenta.hpd.ND.top.rds')

###############
Armenta.hpd.Lambda.top <- list()
Armenta.Lambda.pgls.summary <- na.omit(Armenta.Lambda.pgls.summary)
for(x in unique(Armenta.Lambda.pgls.summary$Parameter)) {
Armenta.hpd.Lambda.top[[x]] = hdi(Armenta.Lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.Lambda.top <- bind_rows(Armenta.hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) 
saveRDS(Armenta.hpd.Lambda.top, 'data/Armenta.hpd.Lambda.top.rds')

################
Armenta.hpd.Mu.top <- list()
Armenta.Mu.pgls.summary <- na.omit(Armenta.Mu.pgls.summary)
for(x in unique(Armenta.Mu.pgls.summary$Parameter)) {
Armenta.hpd.Mu.top[[x]] = hdi(Armenta.Mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.Mu.top <- bind_rows(Armenta.hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Armenta.hpd.Mu.top, 'data/Armenta.hpd.Mu.top.rds')

Armenta.hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta DR HPD Interval")
Armenta DR HPD Interval
  Colour.discriminability log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.0178 -0.0149 -0.000138 0.0107 -0.0425 -4.22e-06
Upper 0.0349 0.0315 0.000161 0.0347 0.00788 1.45e-05
Armenta.hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta ND HPD Interval")
Armenta ND HPD Interval
  Colour.discriminability log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.00176 -0.000596 -4.94e-06 0.000502 -0.0023 -6.87e-08
Upper 0.00173 0.00133 1.38e-05 0.00184 0.000554 6.84e-07
Armenta.hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta BAMM Speciation HPD Interval")
Armenta BAMM Speciation HPD Interval
  Colour.discriminability log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.0127 -0.00712 -4.4e-05 -0.00289 -0.0173 -4.05e-06
Upper 0.011 0.00835 5.99e-05 0.00674 0.00503 3.7e-06
Armenta.hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta BAMM Extinction HPD Interval")
Armenta BAMM Extinction HPD Interval
  Colour.discriminability log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.0282 -0.019 -0.000119 -0.00579 -0.0377 -1.05e-05
Upper 0.0172 0.0206 0.000136 0.0155 0.00783 1.37e-05

Analysis using male-biased measure of sexual selection

Sexual dichromatism is expected to be a good measure of sexual selection strength in birds. However, the relationship between sexual dichromatism and male-biased measures of sexual selection (social mating system, sexual size dimorphism and paternal care) is expected to be relatively noisy given the precision of the measurement used (Dale et al. 2015). Here we use a dataset of sexual selection for 2,465 species of the 5,812 species with sexual dichromatism scores from RGB measures. This male-biased sexual selection score is based on three components, combined in a phylogenetic PCA (ppca) with the following loadings:

  • sexual size dimorphism: 0.37
  • social polygyny: 0.57
  • paternal care: -0.57

As such, species with high values for this score are expected to have high male-biased sexual selection (high dimorphism, high polygyny and low paternal care). PGLS models using this dataset were run using the same process as above.

SS.subset <- plumage.scores %>% filter(Sexual_selection_ppca != "NA")
SS.subset <- inner_join(SS.subset, restricted.data %>% dplyr::select(TipLabel, binomial, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP, MCC.DR, MCC.ND), by = "TipLabel")
SS.subset <- inner_join(SS.subset, MCC.BAMM.df, by = "TipLabel")

SS.subset %>% ggplot(aes(x = SDi, y = Sexual_selection_ppca))+
  geom_point()+
  geom_smooth(method = "lm")+
  theme_minimal()

Figure S12: The relationship between absolute sexual dichromatism and male-baised sexual selection. Although positively correlated, the distribution is very noisy with some species having high sexual dichromatism and low male-biased sexual selection and vice versa.

data_frame(
R2 = summary(lm(SDi ~ Sexual_selection_ppca,
   data = SS.subset))$r.squared,
r = cor(SS.subset$SDi, SS.subset$Sexual_selection_ppca)) %>% pander(digits = 2)
R2 r
0.12 0.34
#Prune tree
pruned.MCC.Subset.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(SS.subset$TipLabel, MCC.passerine$tip.label)])

saveRDS(pruned.MCC.Subset.tree, 'data/pruned.MCC.Subset.tree.rds')

#Set rownames to match tree
rownames(SS.subset) <- SS.subset$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.Subset.global <- gls(MCC.DR ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = FALSE), 
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.ND.Subset.global <- gls(MCC.ND ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), #lambda = 1
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Lambda.Subset.global <- gls(mean.lambda ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variatioColour.discriminability
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), 
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Extinction.Subset.global <- gls(mean.mu ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                weights = ~ sqrt(var.mu),
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), 
                data = SS.subset, 
                method = "REML")
#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("SS.subset"), envir=environment())
clusterExport(myclust, c("pruned.MCC.Subset.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

Subset.dredged.ND.model <- pdredge(MCC.ND.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.ND.model, 'data/Subset.dredged.ND.model.rds')

Subset.dredged.DR.model <- pdredge(MCC.DR.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.DR.model, 'data/Subset.dredged.DR.model.rds')

Subset.dredged.spec.model <- pdredge(MCC.Lambda.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.spec.model, 'data/Subset.dredged.spec.model.rds')

Subset.dredged.extinct.model <- pdredge(MCC.Extinction.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.spec.model, 'data/Subset.dredged.extinct.model.rds')

Table S13: All top models (\(\delta AICc > 4\)) using male-biased sexual selection measures do not contain interactions between sexual selection and the environmental variables. Thus interaction terms were not included in further analysis.

Subset.dredged.ND.model <- readRDS("data/Subset.dredged.ND.model.rds")
Subset.dredged.DR.model <- readRDS("data/Subset.dredged.DR.model.rds")
Subset.dredged.spec.model <- readRDS("data/Subset.dredged.spec.model.rds")
Subset.dredged.extinct.model <- readRDS("data/Subset.dredged.extinct.model.rds")


kable(Subset.dredged.ND.model, "html", caption = "ND Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
ND Male-bias Sexual Selection Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 Sexual_selection_ppca bioclim4:Sexual_selection_ppca log(range.size.m2):Sexual_selection_ppca NPP:Sexual_selection_ppca PC1.LIG:Sexual_selection_ppca residuals.PC1:Sexual_selection_ppca df logLik AICc delta weight
0 0.1484737 -5e-07 -0.0001797 0 0.0001070 0.0000794 0.0004383 NA NA NA NA NA 8 5451.548 -10887.04 0.00000 9.994583e-01
16 0.1482203 -5e-07 -0.0001719 0 0.0000927 0.0000649 0.0003788 NA NA NA NA -0.0001168 9 5444.065 -10870.06 16.98151 2.052470e-04
8 0.1479417 -5e-07 -0.0001568 0 0.0000773 0.0001053 0.0004931 NA NA NA 0.0001281 NA 9 5444.027 -10869.98 17.05724 1.976207e-04
2 0.1483180 -5e-07 -0.0001762 0 0.0001013 0.0000840 0.0017983 NA -4.85e-05 NA NA NA 9 5443.645 -10869.22 17.82173 1.348423e-04
1 0.1478705 -9e-07 -0.0001559 0 0.0000735 0.0000877 0.0012070 -2.3e-06 NA NA NA NA 9 5440.066 -10862.06 24.98007 3.761940e-06
4 0.1484151 -6e-07 -0.0001727 0 0.0001021 0.0000821 0.0008191 NA NA 0e+00 NA NA 9 5436.296 -10854.52 32.51857 8.678530e-08
24 0.1477915 -5e-07 -0.0001526 0 0.0000688 0.0000930 0.0004461 NA NA NA 0.0001198 -0.0000851 10 5436.459 -10852.83 34.21073 3.723916e-08
18 0.1479736 -6e-07 -0.0001660 0 0.0000830 0.0000693 0.0022306 NA -6.64e-05 NA NA -0.0001324 10 5436.227 -10852.37 34.67304 2.955353e-08
10 0.1477977 -5e-07 -0.0001536 0 0.0000721 0.0001095 0.0017750 NA -4.58e-05 NA 0.0001275 NA 10 5436.117 -10852.14 34.89407 2.646140e-08
17 0.1478029 -9e-07 -0.0001541 0 0.0000697 0.0000811 0.0011337 -2.2e-06 NA NA NA -0.0000487 10 5432.460 -10844.83 42.20729 6.832316e-10
9 0.1474756 -8e-07 -0.0001387 0 0.0000514 0.0001091 0.0011730 -2.1e-06 NA NA 0.0001104 NA 10 5432.374 -10844.66 42.38012 6.266710e-10
3 0.1478853 -9e-07 -0.0001561 0 0.0000739 0.0000871 0.0010258 -2.3e-06 6.80e-06 NA NA NA 10 5432.151 -10844.21 42.82649 5.013149e-10
20 0.1482027 -6e-07 -0.0001671 0 0.0000904 0.0000688 0.0006980 NA NA 0e+00 NA -0.0001027 10 5428.778 -10837.47 49.57233 1.718978e-11
12 0.1478903 -6e-07 -0.0001503 0 0.0000728 0.0001075 0.0008573 NA NA 0e+00 0.0001270 NA 10 5428.763 -10837.44 49.60122 1.694322e-11
6 0.1482206 -6e-07 -0.0001678 0 0.0000948 0.0000879 0.0025058 NA -5.89e-05 0e+00 NA NA 10 5428.428 -10836.77 50.27097 1.212167e-11
26 0.1475788 -6e-07 -0.0001477 0 0.0000606 0.0000964 0.0021023 NA -5.94e-05 NA 0.0001176 -0.0000996 11 5428.600 -10835.09 51.94627 5.245357e-12
5 0.1477875 -1e-06 -0.0001473 0 0.0000670 0.0000909 0.0016614 -2.4e-06 NA -1e-07 NA NA 10 5424.858 -10829.63 57.41192 3.411341e-13
25 0.1474471 -8e-07 -0.0001380 0 0.0000498 0.0001053 0.0011354 -2.0e-06 NA NA 0.0001086 -0.0000253 11 5424.753 -10827.40 59.64011 1.119641e-13
19 0.1477901 -9e-07 -0.0001539 0 0.0000692 0.0000812 0.0012557 -2.1e-06 -4.70e-06 NA NA -0.0000505 11 5424.592 -10827.08 59.96212 9.531385e-14
11 0.1474825 -8e-07 -0.0001388 0 0.0000516 0.0001088 0.0010911 -2.1e-06 3.10e-06 NA 0.0001104 NA 11 5424.458 -10826.81 60.22923 8.339781e-14
28 0.1477721 -6e-07 -0.0001476 0 0.0000663 0.0000971 0.0007720 NA NA 0e+00 0.0001202 -0.0000707 11 5421.175 -10820.24 66.79552 3.128260e-15
22 0.1479291 -7e-07 -0.0001601 0 0.0000795 0.0000741 0.0027738 NA -7.31e-05 0e+00 NA -0.0001183 11 5420.967 -10819.83 67.21118 2.541220e-15
14 0.1477101 -6e-07 -0.0001458 0 0.0000661 0.0001129 0.0024524 NA -5.57e-05 0e+00 0.0001261 NA 11 5420.886 -10819.66 67.37372 2.342865e-15
21 0.1477550 -1e-06 -0.0001467 0 0.0000652 0.0000872 0.0016020 -2.3e-06 NA 0e+00 NA -0.0000261 11 5417.250 -10812.39 74.64591 6.174649e-17
13 0.1474029 -1e-06 -0.0001307 0 0.0000456 0.0001118 0.0016079 -2.1e-06 NA 0e+00 0.0001086 NA 11 5417.150 -10812.19 74.84554 5.588074e-17
7 0.1477797 -1e-06 -0.0001472 0 0.0000668 0.0000912 0.0017551 -2.4e-06 -3.40e-06 -1e-07 NA NA 11 5416.949 -10811.79 75.24828 4.568871e-17
27 0.1474393 -8e-07 -0.0001379 0 0.0000496 0.0001054 0.0012104 -2.0e-06 -2.90e-06 NA 0.0001086 -0.0000265 12 5416.884 -10809.64 77.39699 1.560351e-17
30 0.1475332 -6e-07 -0.0001417 0 0.0000570 0.0001013 0.0026495 NA -6.62e-05 0e+00 0.0001178 -0.0000854 12 5413.341 -10802.56 84.48247 4.514702e-19
29 0.1474000 -1e-06 -0.0001307 0 0.0000455 0.0001113 0.0016015 -2.1e-06 NA 0e+00 0.0001084 -0.0000029 12 5409.541 -10794.95 92.08348 1.009462e-20
23 0.1477281 -1e-06 -0.0001462 0 0.0000643 0.0000875 0.0018597 -2.2e-06 -9.80e-06 0e+00 NA -0.0000298 12 5409.385 -10794.64 92.39564 8.635863e-21
15 0.1473872 -1e-06 -0.0001304 0 0.0000451 0.0001124 0.0017903 -2.1e-06 -6.70e-06 0e+00 0.0001087 NA 12 5409.242 -10794.36 92.68146 7.485855e-21
31 0.1473784 -1e-06 -0.0001303 0 0.0000447 0.0001116 0.0018109 -2.1e-06 -7.90e-06 0e+00 0.0001083 -0.0000059 13 5401.675 -10777.20 109.83606 1.409829e-24
kable(Subset.dredged.DR.model, "html", caption = "DR Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
DR Male-bias Sexual Selection Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 Sexual_selection_ppca bioclim4:Sexual_selection_ppca log(range.size.m2):Sexual_selection_ppca NPP:Sexual_selection_ppca PC1.LIG:Sexual_selection_ppca residuals.PC1:Sexual_selection_ppca df logLik AICc delta weight
0 -2.644562 2.44e-05 -0.0112977 -7.0e-07 0.0068773 0.0004563 0.0388704 NA NA NA NA NA 9 -1717.439 3452.952 0.00000 9.851176e-01
16 -2.639831 2.52e-05 -0.0115554 -7.0e-07 0.0066501 0.0000102 0.0381207 NA NA NA NA -0.0043823 10 -1721.449 3462.987 10.03539 6.521251e-03
2 -2.651623 2.26e-05 -0.0111180 -7.0e-07 0.0068694 0.0003268 0.1031437 NA -0.0023145 NA NA NA 10 -1721.904 3463.898 10.94632 4.135475e-03
8 -2.643876 2.47e-05 -0.0113594 -8.0e-07 0.0067897 0.0002975 0.0391199 NA NA NA -0.0017632 NA 10 -1721.917 3463.924 10.97232 4.082050e-03
1 -2.644957 1.72e-05 -0.0112402 -8.0e-07 0.0068984 0.0003008 0.0537607 -4.40e-05 NA NA NA NA 10 -1726.051 3472.191 19.23922 6.542518e-05
18 -2.648516 2.30e-05 -0.0113543 -7.0e-07 0.0066067 -0.0002256 0.1226913 NA -0.0030498 NA NA -0.0050265 11 -1725.773 3473.653 20.70155 3.149221e-05
24 -2.639485 2.54e-05 -0.0115900 -7.0e-07 0.0065973 -0.0000848 0.0383373 NA NA NA -0.0012535 -0.0041982 11 -1725.962 3474.032 21.08051 2.605628e-05
10 -2.650654 2.30e-05 -0.0111829 -7.0e-07 0.0067902 0.0001885 0.1005406 NA -0.0022125 NA -0.0015977 NA 11 -1726.396 3474.899 21.94686 1.689613e-05
4 -2.639330 2.36e-05 -0.0112937 -1.1e-06 0.0067321 0.0005649 0.0550582 NA NA -1.6e-06 NA NA 10 -1729.240 3478.570 25.61774 2.695667e-06
17 -2.641385 1.95e-05 -0.0114418 -8.0e-07 0.0067290 0.0000141 0.0497124 -3.37e-05 NA NA NA -0.0031712 11 -1730.243 3482.593 29.64087 3.606241e-07
9 -2.644482 1.77e-05 -0.0112768 -8.0e-07 0.0068536 0.0002276 0.0532047 -4.20e-05 NA NA -0.0008802 NA 11 -1730.572 3483.251 30.29951 2.594373e-07
3 -2.648691 1.70e-05 -0.0111501 -7.0e-07 0.0068919 0.0002443 0.0875544 -3.97e-05 -0.0012698 NA NA NA 11 -1730.601 3483.309 30.35724 2.520557e-07
26 -2.647952 2.32e-05 -0.0113880 -7.0e-07 0.0065680 -0.0002916 0.1205759 NA -0.0029674 NA -0.0009463 -0.0048704 12 -1730.300 3484.728 31.77614 1.239900e-07
20 -2.635734 2.45e-05 -0.0115337 -1.0e-06 0.0065447 0.0001349 0.0518435 NA NA -1.4e-06 NA -0.0040546 11 -1733.313 3488.734 35.78248 1.672712e-08
6 -2.647032 2.13e-05 -0.0110757 -1.1e-06 0.0066999 0.0004272 0.1348340 NA -0.0027804 -1.9e-06 NA NA 11 -1733.618 3489.344 36.39193 1.233336e-08
12 -2.638528 2.39e-05 -0.0113579 -1.1e-06 0.0066376 0.0004005 0.0556204 NA NA -1.7e-06 -0.0018504 NA 11 -1733.708 3489.524 36.57215 1.127059e-08
19 -2.647369 1.96e-05 -0.0113287 -7.0e-07 0.0066708 -0.0001664 0.1099449 -2.28e-05 -0.0023064 NA NA -0.0040490 12 -1734.642 3493.411 40.45959 1.613617e-09
25 -2.641013 1.98e-05 -0.0114694 -8.0e-07 0.0066956 -0.0000410 0.0493348 -3.22e-05 NA NA -0.0007088 -0.0031213 12 -1734.769 3493.664 40.71256 1.421894e-09
11 -2.648207 1.74e-05 -0.0111870 -8.0e-07 0.0068471 0.0001713 0.0869569 -3.76e-05 -0.0012682 NA -0.0008789 NA 12 -1735.121 3494.370 41.41788 9.993313e-10
5 -2.639820 1.65e-05 -0.0112370 -1.1e-06 0.0067562 0.0004086 0.0694696 -4.36e-05 NA -1.6e-06 NA NA 11 -1737.861 3497.830 44.87795 1.771599e-10
22 -2.644671 2.18e-05 -0.0113032 -1.1e-06 0.0064753 -0.0001033 0.1491840 NA -0.0034111 -1.7e-06 NA -0.0047084 12 -1737.558 3499.243 46.29127 8.739096e-11
28 -2.635275 2.46e-05 -0.0115702 -1.1e-06 0.0064843 0.0000345 0.0524329 NA NA -1.4e-06 -0.0013712 -0.0038442 12 -1737.817 3499.761 46.80906 6.745731e-11
14 -2.645993 2.16e-05 -0.0111425 -1.2e-06 0.0066161 0.0002845 0.1324246 NA -0.0026787 -1.9e-06 -0.0016609 NA 12 -1738.103 3500.333 47.38132 5.067171e-11
27 -2.646963 2.00e-05 -0.0113557 -7.0e-07 0.0066402 -0.0002164 0.1092438 -2.15e-05 -0.0022928 NA -0.0006534 -0.0039979 13 -1739.169 3504.487 51.53478 6.351151e-12
21 -2.637174 1.86e-05 -0.0114163 -1.1e-06 0.0066217 0.0001435 0.0642405 -3.45e-05 NA -1.4e-06 NA -0.0028007 12 -1742.094 3508.315 55.36277 9.367329e-13
7 -2.644544 1.60e-05 -0.0111093 -1.1e-06 0.0067323 0.0003423 0.1180262 -3.75e-05 -0.0017625 -1.8e-06 NA NA 12 -1742.353 3508.834 55.88233 7.224286e-13
13 -2.639263 1.70e-05 -0.0112769 -1.1e-06 0.0067049 0.0003287 0.0690245 -4.13e-05 NA -1.6e-06 -0.0009783 NA 12 -1742.376 3508.880 55.92803 7.061065e-13
30 -2.644024 2.20e-05 -0.0113390 -1.1e-06 0.0064309 -0.0001745 0.1471552 NA -0.0033240 -1.7e-06 -0.0010472 -0.0045312 13 -1742.079 3510.306 57.35440 3.460502e-13
23 -2.643615 1.86e-05 -0.0112795 -1.1e-06 0.0065390 -0.0000483 0.1365076 -2.20e-05 -0.0026899 -1.6e-06 NA -0.0037721 13 -1746.434 3519.017 66.06483 4.443157e-15
29 -2.636728 1.90e-05 -0.0114468 -1.1e-06 0.0065818 0.0000821 0.0639906 -3.28e-05 NA -1.4e-06 -0.0008185 -0.0027377 13 -1746.614 3519.376 66.42447 3.711904e-15
15 -2.643986 1.65e-05 -0.0111494 -1.1e-06 0.0066807 0.0002619 0.1176664 -3.52e-05 -0.0017657 -1.8e-06 -0.0009846 NA 13 -1746.868 3519.884 66.93197 2.880021e-15
31 -2.643141 1.90e-05 -0.0113095 -1.1e-06 0.0065019 -0.0001051 0.1359381 -2.04e-05 -0.0026772 -1.6e-06 -0.0007673 -0.0037084 14 -1750.955 3530.082 77.13061 1.757073e-17
kable(Subset.dredged.spec.model, "html", caption = "BAMM Speciation Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM Speciation Male-bias Sexual Selection Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 Sexual_selection_ppca bioclim4:Sexual_selection_ppca log(range.size.m2):Sexual_selection_ppca NPP:Sexual_selection_ppca PC1.LIG:Sexual_selection_ppca residuals.PC1:Sexual_selection_ppca df logLik AICc delta weight
0 -2.291078 -7.6e-06 0.0005462 -1e-07 -0.0005682 0.0004004 0.0009423 NA NA NA NA NA 8 2239.622 -4463.186 0.00000 9.982422e-01
16 -2.291662 -7.7e-06 0.0005642 -1e-07 -0.0005920 0.0003587 0.0008223 NA NA NA NA -0.0002634 9 2233.362 -4448.650 14.53595 6.962976e-04
8 -2.292324 -7.7e-06 0.0005956 -1e-07 -0.0006033 0.0004432 0.0009639 NA NA NA 0.0003153 NA 9 2233.112 -4448.150 15.03549 5.424008e-04
2 -2.290956 -7.5e-06 0.0005410 -1e-07 -0.0005536 0.0003878 -0.0024879 NA 0.0001233 NA NA NA 9 2233.053 -4448.033 15.15309 5.114272e-04
1 -2.291377 -7.8e-06 0.0005595 -1e-07 -0.0005821 0.0004069 0.0013435 -1.2e-06 NA NA NA NA 9 2228.665 -4439.256 23.93013 6.351456e-06
24 -2.292789 -7.8e-06 0.0006095 -1e-07 -0.0006228 0.0004052 0.0008583 NA NA NA 0.0003040 -0.0002302 10 2226.838 -4433.586 29.60044 3.728916e-07
18 -2.291526 -7.6e-06 0.0005589 -1e-07 -0.0005793 0.0003527 -0.0016890 NA 0.0000906 NA NA -0.0002427 10 2226.796 -4433.503 29.68333 3.577523e-07
4 -2.291631 -7.9e-06 0.0005731 -1e-07 -0.0005919 0.0004057 0.0023363 NA NA -2e-07 NA NA 9 2225.686 -4433.299 29.88725 3.230739e-07
10 -2.292200 -7.6e-06 0.0005903 -1e-07 -0.0005890 0.0004308 -0.0023426 NA 0.0001188 NA 0.0003140 NA 10 2226.541 -4432.992 30.19348 2.772079e-07
17 -2.291726 -7.8e-06 0.0005675 -1e-07 -0.0005951 0.0003631 0.0009581 -4.0e-07 NA NA NA -0.0002489 10 2222.455 -4424.820 38.36631 4.656914e-09
3 -2.291372 -7.9e-06 0.0005595 -1e-07 -0.0005694 0.0003930 -0.0032228 -1.9e-06 0.0001721 NA NA NA 10 2222.168 -4424.247 38.93918 3.497033e-09
9 -2.292553 -7.9e-06 0.0006058 -1e-07 -0.0006143 0.0004480 0.0012916 -1.0e-06 NA NA 0.0003118 NA 10 2222.150 -4424.211 38.97537 3.434328e-09
20 -2.292050 -7.9e-06 0.0005850 -1e-07 -0.0006088 0.0003712 0.0020972 NA NA -1e-07 NA -0.0002145 10 2219.414 -4418.739 44.44697 2.226889e-10
26 -2.292653 -7.7e-06 0.0006042 -1e-07 -0.0006101 0.0003992 -0.0016579 NA 0.0000908 NA 0.0003040 -0.0002094 11 2220.272 -4418.437 44.74913 1.914638e-10
12 -2.292849 -8.0e-06 0.0006212 -1e-07 -0.0006260 0.0004479 0.0023254 NA NA -1e-07 0.0003116 NA 10 2219.171 -4418.252 44.93438 1.745257e-10
6 -2.291516 -7.8e-06 0.0005680 -1e-07 -0.0005801 0.0003963 -0.0002245 NA 0.0000897 -1e-07 NA NA 10 2219.114 -4418.138 45.04768 1.649133e-10
19 -2.291648 -7.8e-06 0.0005658 -1e-07 -0.0005832 0.0003623 -0.0022522 -1.1e-06 0.0001241 NA NA -0.0001960 11 2215.990 -4409.872 53.31398 2.643940e-12
25 -2.292831 -7.8e-06 0.0006117 -1e-07 -0.0006249 0.0004082 0.0009520 -3.0e-07 NA NA 0.0003035 -0.0002202 11 2215.930 -4409.753 53.43304 2.491138e-12
5 -2.291991 -8.2e-06 0.0005892 -1e-07 -0.0006086 0.0004135 0.0028374 -1.4e-06 NA -2e-07 NA NA 10 2214.736 -4409.382 53.80420 2.069196e-12
11 -2.292533 -8.0e-06 0.0006053 -1e-07 -0.0006019 0.0004344 -0.0029783 -1.7e-06 0.0001610 NA 0.0003078 NA 11 2215.649 -4409.190 53.99621 1.879781e-12
28 -2.293170 -8.0e-06 0.0006300 -1e-07 -0.0006394 0.0004175 0.0021230 NA NA -1e-07 0.0003030 -0.0001818 11 2212.889 -4403.670 59.51566 1.190073e-13
22 -2.291940 -7.9e-06 0.0005806 -1e-07 -0.0005991 0.0003665 0.0002414 NA 0.0000655 -1e-07 NA -0.0002010 11 2212.849 -4403.590 59.59632 1.143037e-13
14 -2.292735 -7.9e-06 0.0006163 -1e-07 -0.0006146 0.0004387 -0.0001310 NA 0.0000860 -1e-07 0.0003108 NA 11 2212.598 -4403.088 60.09824 8.893407e-14
21 -2.292190 -8.1e-06 0.0005922 -1e-07 -0.0006157 0.0003804 0.0024089 -8.0e-07 NA -1e-07 NA -0.0001840 11 2208.520 -4394.932 68.25439 1.506547e-15
27 -2.292752 -7.9e-06 0.0006100 -1e-07 -0.0006133 0.0004073 -0.0021442 -9.0e-07 0.0001197 NA 0.0003025 -0.0001693 12 2209.464 -4394.801 68.38491 1.411373e-15
7 -2.291949 -8.2e-06 0.0005874 -1e-07 -0.0005966 0.0004017 -0.0009573 -2.0e-06 0.0001396 -1e-07 NA NA 11 2208.231 -4394.355 68.83076 1.129342e-15
13 -2.293134 -8.2e-06 0.0006341 -1e-07 -0.0006396 0.0004538 0.0027470 -1.2e-06 NA -2e-07 0.0003072 NA 11 2208.215 -4394.323 68.86294 1.111321e-15
30 -2.293059 -8.0e-06 0.0006256 -1e-07 -0.0006297 0.0004128 0.0002565 NA 0.0000659 -1e-07 0.0003030 -0.0001682 12 2206.323 -4388.520 74.66625 6.104709e-17
23 -2.292110 -8.1e-06 0.0005902 -1e-07 -0.0006048 0.0003793 -0.0004064 -1.4e-06 0.0001073 -1e-07 NA -0.0001400 12 2202.052 -4379.976 83.20972 8.520698e-19
29 -2.293283 -8.1e-06 0.0006359 -1e-07 -0.0006450 0.0004250 0.0023847 -7.0e-07 NA -1e-07 0.0003018 -0.0001563 12 2201.993 -4379.859 83.32717 8.034705e-19
15 -2.293085 -8.3e-06 0.0006320 -1e-07 -0.0006282 0.0004425 -0.0007674 -1.7e-06 0.0001293 -1e-07 0.0003043 NA 12 2201.707 -4379.288 83.89838 6.038562e-19
31 -2.293203 -8.2e-06 0.0006338 -1e-07 -0.0006345 0.0004238 -0.0003206 -1.2e-06 0.0001031 -1e-07 0.0003009 -0.0001140 13 2195.524 -4364.900 98.28616 4.535936e-22
kable(Subset.dredged.extinct.model, "html", caption = "BAMM Extinction Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM Extinction Male-bias Sexual Selection Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 Sexual_selection_ppca bioclim4:Sexual_selection_ppca log(range.size.m2):Sexual_selection_ppca NPP:Sexual_selection_ppca PC1.LIG:Sexual_selection_ppca residuals.PC1:Sexual_selection_ppca df logLik AICc delta weight
0 -2.291078 -7.6e-06 0.0005462 -1e-07 -0.0005682 0.0004004 0.0009423 NA NA NA NA NA 8 2239.622 -4463.186 0.00000 9.982422e-01
16 -2.291662 -7.7e-06 0.0005642 -1e-07 -0.0005920 0.0003587 0.0008223 NA NA NA NA -0.0002634 9 2233.362 -4448.650 14.53595 6.962976e-04
8 -2.292324 -7.7e-06 0.0005956 -1e-07 -0.0006033 0.0004432 0.0009639 NA NA NA 0.0003153 NA 9 2233.112 -4448.150 15.03549 5.424008e-04
2 -2.290956 -7.5e-06 0.0005410 -1e-07 -0.0005536 0.0003878 -0.0024879 NA 0.0001233 NA NA NA 9 2233.053 -4448.033 15.15309 5.114272e-04
1 -2.291377 -7.8e-06 0.0005595 -1e-07 -0.0005821 0.0004069 0.0013435 -1.2e-06 NA NA NA NA 9 2228.665 -4439.256 23.93013 6.351456e-06
24 -2.292789 -7.8e-06 0.0006095 -1e-07 -0.0006228 0.0004052 0.0008583 NA NA NA 0.0003040 -0.0002302 10 2226.838 -4433.586 29.60044 3.728916e-07
18 -2.291526 -7.6e-06 0.0005589 -1e-07 -0.0005793 0.0003527 -0.0016890 NA 0.0000906 NA NA -0.0002427 10 2226.796 -4433.503 29.68333 3.577523e-07
4 -2.291631 -7.9e-06 0.0005731 -1e-07 -0.0005919 0.0004057 0.0023363 NA NA -2e-07 NA NA 9 2225.686 -4433.299 29.88725 3.230739e-07
10 -2.292200 -7.6e-06 0.0005903 -1e-07 -0.0005890 0.0004308 -0.0023426 NA 0.0001188 NA 0.0003140 NA 10 2226.541 -4432.992 30.19348 2.772079e-07
17 -2.291726 -7.8e-06 0.0005675 -1e-07 -0.0005951 0.0003631 0.0009581 -4.0e-07 NA NA NA -0.0002489 10 2222.455 -4424.820 38.36631 4.656914e-09
3 -2.291372 -7.9e-06 0.0005595 -1e-07 -0.0005694 0.0003930 -0.0032228 -1.9e-06 0.0001721 NA NA NA 10 2222.168 -4424.247 38.93918 3.497033e-09
9 -2.292553 -7.9e-06 0.0006058 -1e-07 -0.0006143 0.0004480 0.0012916 -1.0e-06 NA NA 0.0003118 NA 10 2222.150 -4424.211 38.97537 3.434328e-09
20 -2.292050 -7.9e-06 0.0005850 -1e-07 -0.0006088 0.0003712 0.0020972 NA NA -1e-07 NA -0.0002145 10 2219.414 -4418.739 44.44697 2.226889e-10
26 -2.292653 -7.7e-06 0.0006042 -1e-07 -0.0006101 0.0003992 -0.0016579 NA 0.0000908 NA 0.0003040 -0.0002094 11 2220.272 -4418.437 44.74913 1.914638e-10
12 -2.292849 -8.0e-06 0.0006212 -1e-07 -0.0006260 0.0004479 0.0023254 NA NA -1e-07 0.0003116 NA 10 2219.171 -4418.252 44.93438 1.745257e-10
6 -2.291516 -7.8e-06 0.0005680 -1e-07 -0.0005801 0.0003963 -0.0002245 NA 0.0000897 -1e-07 NA NA 10 2219.114 -4418.138 45.04768 1.649133e-10
19 -2.291648 -7.8e-06 0.0005658 -1e-07 -0.0005832 0.0003623 -0.0022522 -1.1e-06 0.0001241 NA NA -0.0001960 11 2215.990 -4409.872 53.31398 2.643940e-12
25 -2.292831 -7.8e-06 0.0006117 -1e-07 -0.0006249 0.0004082 0.0009520 -3.0e-07 NA NA 0.0003035 -0.0002202 11 2215.930 -4409.753 53.43304 2.491138e-12
5 -2.291991 -8.2e-06 0.0005892 -1e-07 -0.0006086 0.0004135 0.0028374 -1.4e-06 NA -2e-07 NA NA 10 2214.736 -4409.382 53.80420 2.069196e-12
11 -2.292533 -8.0e-06 0.0006053 -1e-07 -0.0006019 0.0004344 -0.0029783 -1.7e-06 0.0001610 NA 0.0003078 NA 11 2215.649 -4409.190 53.99621 1.879781e-12
28 -2.293170 -8.0e-06 0.0006300 -1e-07 -0.0006394 0.0004175 0.0021230 NA NA -1e-07 0.0003030 -0.0001818 11 2212.889 -4403.670 59.51566 1.190073e-13
22 -2.291940 -7.9e-06 0.0005806 -1e-07 -0.0005991 0.0003665 0.0002414 NA 0.0000655 -1e-07 NA -0.0002010 11 2212.849 -4403.590 59.59632 1.143037e-13
14 -2.292735 -7.9e-06 0.0006163 -1e-07 -0.0006146 0.0004387 -0.0001310 NA 0.0000860 -1e-07 0.0003108 NA 11 2212.598 -4403.088 60.09824 8.893407e-14
21 -2.292190 -8.1e-06 0.0005922 -1e-07 -0.0006157 0.0003804 0.0024089 -8.0e-07 NA -1e-07 NA -0.0001840 11 2208.520 -4394.932 68.25439 1.506547e-15
27 -2.292752 -7.9e-06 0.0006100 -1e-07 -0.0006133 0.0004073 -0.0021442 -9.0e-07 0.0001197 NA 0.0003025 -0.0001693 12 2209.464 -4394.801 68.38491 1.411373e-15
7 -2.291949 -8.2e-06 0.0005874 -1e-07 -0.0005966 0.0004017 -0.0009573 -2.0e-06 0.0001396 -1e-07 NA NA 11 2208.231 -4394.355 68.83076 1.129342e-15
13 -2.293134 -8.2e-06 0.0006341 -1e-07 -0.0006396 0.0004538 0.0027470 -1.2e-06 NA -2e-07 0.0003072 NA 11 2208.215 -4394.323 68.86294 1.111321e-15
30 -2.293059 -8.0e-06 0.0006256 -1e-07 -0.0006297 0.0004128 0.0002565 NA 0.0000659 -1e-07 0.0003030 -0.0001682 12 2206.323 -4388.520 74.66625 6.104709e-17
23 -2.292110 -8.1e-06 0.0005902 -1e-07 -0.0006048 0.0003793 -0.0004064 -1.4e-06 0.0001073 -1e-07 NA -0.0001400 12 2202.052 -4379.976 83.20972 8.520698e-19
29 -2.293283 -8.1e-06 0.0006359 -1e-07 -0.0006450 0.0004250 0.0023847 -7.0e-07 NA -1e-07 0.0003018 -0.0001563 12 2201.993 -4379.859 83.32717 8.034705e-19
15 -2.293085 -8.3e-06 0.0006320 -1e-07 -0.0006282 0.0004425 -0.0007674 -1.7e-06 0.0001293 -1e-07 0.0003043 NA 12 2201.707 -4379.288 83.89838 6.038562e-19
31 -2.293203 -8.2e-06 0.0006338 -1e-07 -0.0006345 0.0004238 -0.0003206 -1.2e-06 0.0001031 -1e-07 0.0003009 -0.0001140 13 2195.524 -4364.900 98.28616 4.535936e-22
#Run model for DR
Subset.MCC.DR.top <- gls(MCC.DR ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = FALSE), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.DR.top, 'data/Subset.MCC.DR.top.rds')

#Run model for ND
Subset.MCC.ND.top <- gls(MCC.ND ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), #lambda = 1
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.ND.top, 'data/Subset.MCC.ND.top.rds')

## No longer needed 

# #Run the 100 models for DR and ND using the best model:
# Subset.data.noMCC <- SS.subset %>% dplyr::select(TipLabel, Sexual_selection_ppca, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP)
# 
# #Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
# Subset.DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
#   left_join(Subset.data.noMCC, 
#             x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "DR")), 
#             by = "TipLabel")
# })
# 
# #PGLS needs tiplabel as rowname
# Subset.DR.model.data <- lapply(Subset.DR.model.data, function(x) {
#   tibble::column_to_rownames(x, "TipLabel")})
# 
# #Prune the trees
# Subset.pruned.trees<-lapply(passerine.trees, function(x) {
#   drop.tip(x,x$tip.label[-match(SS.subset$TipLabel, x$tip.label)])
# })
# 
# #Use mapply to create a list of PGLS global models
# Subset.DR.pgls.models <- mcmapply(function(x,y) {
#   gls(DR ~ Sexual_selection_ppca 
#          + log(range.size.m2)
#          + bioclim4 #Seasonal variation
#          + residuals.PC1 #Spatial variation
#          + PC1.LIG #Long-term climate variation
#          + NPP,
#     correlation = corPagel(0.9711, phy = y, fixed = TRUE), 
#     data = x, 
#     method = "REML")
# }, x = Subset.DR.model.data, y = Subset.pruned.trees,
# SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
# mc.cores = 8) #Specify core number 
# 
# saveRDS(Subset.DR.pgls.models, "data/Subset.DR.pgls.models.rds")
# 
# #Now for Node Density:
# Subset.ND.model.data <- lapply(nd.list, function(x) { #es.list is a list of ES values calculated earlier
#   left_join(Subset.data.noMCC, 
#             x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
#             by = "TipLabel")
# })
# 
# #PGLS needs tiplabel as rowname
# Subset.ND.model.data <- lapply(Subset.ND.model.data, function(x) {
#   tibble::column_to_rownames(x, "TipLabel")})
# 
# #Use mapply to create a list of PGLS global models
# Subset.ND.pgls.models <- mcmapply(function(x,y) {
# gls(ND ~ Sexual_selection_ppca 
#          + log(range.size.m2)
#          + bioclim4 #Seasonal variation
#          + residuals.PC1 #Spatial variation
#          + PC1.LIG #Long-term climate variation
#          + NPP,
#     corPagel(1, phy = y, fixed = TRUE), 
#     data = x, 
#     method = "REML")
# }, x = Subset.ND.model.data, y = Subset.pruned.trees, 
# SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
# mc.cores = 8) #Specify core number 
# 
# saveRDS(Subset.ND.pgls.models, "data/Subset.ND.pgls.models.rds")

#BAMM Top Models 

Subset.MCC.Lambda.top <- gls(mean.lambda ~ Sexual_selection_ppca
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corBrownian(phy = pruned.MCC.Subset.tree), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.Lambda.top, 'data/Subset.MCC.Lambda.top.rds')

Subset.MCC.Mu.top <- gls(mean.mu ~ Sexual_selection_ppca
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.Subset.tree), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.Mu.top, 'data/Subset.MCC.Mu.top.rds')

#Now for BAMM 100 models

Subset.BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Subset.data.noMCC %>% dplyr::select(TipLabel, Sexual_selection_ppca, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Subset.BAMM.model.data <- lapply(Subset.BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Subset.BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(mean.lambda ~ Sexual_selection_ppca 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.lambda),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Subset.BAMM.model.data, y = Subset.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.BAMM.lambda.pgls.models, "data/Subset.BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
Subset.BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(mean.mu ~ Sexual_selection_ppca
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.mu),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Subset.BAMM.model.data, y = Subset.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.BAMM.mu.pgls.models, "data/Subset.BAMM.mu.pgls.models.rds")

The phylogenetic signal from the DR model is ~0.97. The other models had phylogenetic signals ~ 1, as such they were fixed at 1 to avoid problems of convergence.

Subset.MCC.DR.top <- readRDS('data/Subset.MCC.DR.top.rds')
Subset.MCC.DR.top[["modelStruct"]][["corStruct"]] %>% `names<-`("DR lambda") %>% pander()
  • DR lambda: 0.9711
##########################################################
####### Read in DR models and extract estimates: #########
##########################################################

#We ran DR and ND models on 1000 trees
files.DR.SS <- list.files(path = "/Users/justincally/Dropbox/Runs Spartan/DR_SS/", pattern = "\\.rds$", full.names = TRUE) #1000 models
df <- list()
rds.list <- list()
models <- list()
c.list <- NULL

lapply(files.DR.SS, function(x) {
models <- readRDS(x)
rds.list <- unlist(models, recursive = F, use.names = T)
c.list <<- c(c.list, rds.list[-grep(".log", names(rds.list))])
rm(models)
rm(rds.list)
gc()
})


lapply(names(c.list), function(x) {
  tryCatch({
  model <- c.list[[x]]
  df[[x]] <<- data.frame(model$coefficients,
                              confint(model),
                              coef(summary(model))[,2], #Std.Error
                              coef(summary(model))[,3], #t-val
                              coef(summary(model))[,4], #pval
                              model[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = x) %>% tibble::rownames_to_column()
  gc(verbose = FALSE)
  },
  error = function(e) NULL
  )
})

Subset.DR.pgls.summary <- bind_rows(df) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name"))
#saveRDS(Subset.DR.pgls.summary, "data/Subset.DR.pgls.summary.rds") #Save a simple df with the model and coeff



#read in 1000 trees: 

#We ran DR and ND models on 1000 trees
files.ND.SS <- list.files(path = "/Users/justincally/Dropbox/Runs Spartan/ND_SS/", pattern = "\\.rds$", full.names = TRUE) #1000 models
df <- list()
rds.list <- list()
models <- list()
c.list <- NULL

lapply(files.ND.SS, function(x) {
models <- readRDS(x)
rds.list <- unlist(models, recursive = F, use.names = T)
c.list <<- c(c.list, rds.list[-grep(".log", names(rds.list))])
rm(models)
rm(rds.list)
gc()
})


lapply(names(c.list), function(x) {
  tryCatch({
  model <- c.list[[x]]
  df[[x]] <<- data.frame(model$coefficients,
                              confint(model),
                              coef(summary(model))[,2], #Std.Error
                              coef(summary(model))[,3], #t-val
                              coef(summary(model))[,4], #pval
                              model[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = x) %>% tibble::rownames_to_column()
  gc(verbose = FALSE)
  },
  error = function(e) NULL
  )
})

Subset.ND.pgls.summary <- bind_rows(df) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name"))
#saveRDS(Subset.ND.pgls.summary, "data/Subset.ND.pgls.summary.rds") #Save a simple df with the model and coeff
rm(c.list)
############################
Subset.DR.pgls.summary <- readRDS("data/Subset.DR.pgls.summary.rds")
Subset.ND.pgls.summary <- readRDS("data/Subset.ND.pgls.summary.rds")

#Cols for plot: 
Subset.cols <- brewer.pal(n = 7, name = "Dark2")[-6]

### MCC
Subset.MCC.DR.top <- readRDS('data/Subset.MCC.DR.top.rds')

Subset.MCC.DR.summary <- data.frame(Subset.MCC.DR.top$coefficients,
                              confint(Subset.MCC.DR.top),
                              coef(summary(Subset.MCC.DR.top))[,2], #Std.Error
                              coef(summary(Subset.MCC.DR.top))[,3], #t-val
                              coef(summary(Subset.MCC.DR.top))[,4], #pval
                              Subset.MCC.DR.top[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = "MCC_model") %>% tibble::rownames_to_column()

colnames(Subset.MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `Sexual_selection_ppca` = "Sexual Selection"
                    )

Subset.DR.pgls.summary$Parameter = factor(Subset.DR.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.DR.summary$Parameter = factor(Subset.MCC.DR.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.DR.plot <-Subset.DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.DR.plot <- Subset.DR.plot + geom_errorbarh(data = Subset.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS DR")

####################################
#For ND

Subset.MCC.ND.top <- readRDS('data/Subset.MCC.ND.top.rds')

Subset.MCC.ND.summary <- data.frame(Subset.MCC.ND.top$coefficients,
                              confint(Subset.MCC.ND.top),
                              coef(summary(Subset.MCC.ND.top))[,2], #Std.Error
                              coef(summary(Subset.MCC.ND.top))[,3], #t-val
                              coef(summary(Subset.MCC.ND.top))[,4], #pval
                              Subset.MCC.ND.top[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = "MCC_model") %>% tibble::rownames_to_column()

colnames(Subset.MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name")

Subset.ND.pgls.summary$Parameter = factor(Subset.ND.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.ND.summary$Parameter = factor(Subset.MCC.ND.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.ND.plot <-Subset.ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.ND.plot <- Subset.ND.plot + geom_errorbarh(data = Subset.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS ND")

######################

#For Lambda
Subset.BAMM.lambda.pgls.models <- readRDS("data/Subset.BAMM.lambda.pgls.models.rds")
Subset.MCC.Lambda.top <- readRDS('data/Subset.MCC.Lambda.top.rds')

Subset.Lambda.pgls.summary <- lapply(Subset.BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.Lambda.pgls.summary <- bind_rows(Subset.Lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

#MCC

Subset.MCC.Lambda.summary <- data.frame(Subset.MCC.Lambda.top$coefficients,
                              confint(Subset.MCC.Lambda.top),
                              coef(summary(Subset.MCC.Lambda.top))[,2], #Std.Error
                              coef(summary(Subset.MCC.Lambda.top))[,3], #t-val
                              coef(summary(Subset.MCC.Lambda.top))[,4], #pval
                              Subset.MCC.Lambda.top[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = "MCC_model") %>% tibble::rownames_to_column()

colnames(Subset.MCC.Lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name")


Subset.Lambda.pgls.summary$Parameter = factor(Subset.Lambda.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.Lambda.summary$Parameter = factor(Subset.MCC.Lambda.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.Lambda.pgls.summary.RO <- dcast(Subset.Lambda.pgls.summary %>% filter(Parameter != "(Intercept)"), Estimate ~ Parameter, value.var = "Estimate")
Subset.Lambda.pgls.summary.RO$Estimate <- NULL
Subset.Lambda.pgls.summary.RO <- sapply(Subset.Lambda.pgls.summary.RO, function(x) {
  remove_outliers(x, na.rm = T)})
Subset.Lambda.pgls.summary.RO <-melt(Subset.Lambda.pgls.summary.RO) %>% na.omit()
Subset.Lambda.pgls.summary.RO$Var1 <- NULL
colnames(Subset.Lambda.pgls.summary.RO) <- c("Parameter", "Estimate")


Subset.Lambda.plot <-Subset.Lambda.pgls.summary.RO %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.Lambda.plot <- Subset.Lambda.plot + geom_errorbarh(data = Subset.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS Speciaton")

######################

#For Mu

Subset.BAMM.Mu.pgls.models <- readRDS("data/Subset.BAMM.mu.pgls.models.rds")
Subset.MCC.Mu.top <- readRDS('data/Subset.MCC.Mu.top.rds')

Subset.Mu.pgls.summary <- lapply(Subset.BAMM.Mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})
Subset.Mu.pgls.summary <- bind_rows(Subset.Mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

#MCC

Subset.MCC.Mu.summary <- data.frame(Subset.MCC.Mu.top$coefficients,
                              confint(Subset.MCC.Mu.top),
                              coef(summary(Subset.MCC.Mu.top))[,2], #Std.Error
                              coef(summary(Subset.MCC.Mu.top))[,3], #t-val
                              coef(summary(Subset.MCC.Mu.top))[,4], #pval
                              Subset.MCC.Mu.top[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = "MCC_model") %>% tibble::rownames_to_column()

colnames(Subset.MCC.Mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name")

Subset.Mu.pgls.summary$Parameter = factor(Subset.Mu.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.Mu.summary$Parameter = factor(Subset.MCC.Mu.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))


Subset.Mu.pgls.summary.RO <- dcast(Subset.Mu.pgls.summary %>% filter(Parameter != "(Intercept)"), Estimate ~ Parameter, value.var = "Estimate")
Subset.Mu.pgls.summary.RO$Estimate <- NULL
Subset.Mu.pgls.summary.RO <- sapply(Subset.Mu.pgls.summary.RO, function(x) {
  remove_outliers(x, na.rm = T)})
Subset.Mu.pgls.summary.RO <-melt(Subset.Mu.pgls.summary.RO) %>% na.omit()
Subset.Mu.pgls.summary.RO$Var1 <- NULL
colnames(Subset.Mu.pgls.summary.RO) <- c("Parameter", "Estimate")

Subset.Mu.plot <-Subset.Mu.pgls.summary.RO %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.Mu.plot <- Subset.Mu.plot + geom_errorbarh(data = Subset.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS Extinction")

#########################################
grid.arrange(symmetrise_scale(Subset.DR.plot, "x"),
             symmetrise_scale(Subset.ND.plot, "x"),
             symmetrise_scale(Subset.Lambda.plot, "x"), 
             symmetrise_scale(Subset.Mu.plot, "x"), 
             nrow = 1)

Figure S13: PGLS analyses using measures of male-biased sexual selection (n = 2,465) using three measures of speciation (\(\lambda_{DR}\), \(\lambda_{ND}\), \(\lambda_{BAMM}\)) and one measure of extinction (\(\mu_{BAMM}\)) as response variables. The numerical values for the model estimates using the MCC tree and HPD intervals of estimates from 100 random trees can be found in the ESM. Density curves are based on estimates from 100 trees and the circle below with error bars is the estimate and 95 % CIs from the MCC tree. For this figure we removed outliers from estimates coming from the 100 random trees for BAMM models in order to interpret the MCC 95 % CIs. This figure is part b of Figure 1 within the manuscript.

Table S14: MCC model estimates from the above model are presented here as numberical values with 95 % CIs.

Subset.MCC.DR.summary %>% select(-model_name) %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC DR Estimates from Male-bias SS", split.table = Inf)
MCC DR Estimates from Male-bias SS
Parameter Estimate LCI UCI SE tval pval lambda
Sexual_selection_ppca 0.03887 0.00947 0.06827 0.015 2.591 0.009618 0.9711
log(range.size.m2) -0.0113 -0.02054 -0.002052 0.004718 -2.395 0.0167 0.9711
bioclim4 2.444e-05 -6.115e-05 0.00011 4.367e-05 0.5596 0.5758 0.9711
residuals.PC1 0.0004563 -0.008575 0.009487 0.004608 0.09903 0.9211 0.9711
PC1.LIG 0.006877 -0.001672 0.01543 0.004362 1.577 0.115 0.9711
NPP -7.331e-07 -3.45e-06 1.984e-06 1.386e-06 -0.5289 0.5969 0.9711
Subset.MCC.ND.summary %>% select(-model_name) %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC ND Estimates from Male-bias SS", split.table = Inf)
MCC ND Estimates from Male-bias SS
Parameter Estimate LCI UCI SE tval pval lambda
Sexual_selection_ppca 0.0004383 -0.0004175 0.001294 0.0004367 1.004 0.3156 1
log(range.size.m2) -0.0001797 -0.0004833 0.0001239 0.0001549 -1.16 0.246 1
bioclim4 -4.598e-07 -3.598e-06 2.679e-06 1.601e-06 -0.2872 0.774 1
residuals.PC1 7.944e-05 -0.0002525 0.0004114 0.0001694 0.4691 0.6391 1
PC1.LIG 0.000107 -0.0001581 0.0003721 0.0001352 0.791 0.429 1
NPP -1.722e-08 -1.229e-07 8.841e-08 5.389e-08 -0.3196 0.7493 1
Subset.MCC.Lambda.summary %>% select(-model_name) %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Speciation Estimates from Male-bias SS", split.table = Inf)
MCC BAMM Speciation Estimates from Male-bias SS
Parameter Estimate LCI UCI SE tval pval lambda
Sexual_selection_ppca 0.0009423 -0.002845 0.004729 0.001932 0.4877 0.6258 1
log(range.size.m2) 0.0005462 -0.0006169 0.001709 0.0005934 0.9204 0.3575 1
bioclim4 -7.591e-06 -1.936e-05 4.175e-06 6.003e-06 -1.264 0.2062 1
residuals.PC1 0.0004004 -0.0008464 0.001647 0.0006361 0.6294 0.5291 1
PC1.LIG -0.0005682 -0.001626 0.0004899 0.0005399 -1.053 0.2927 1
NPP -5.972e-08 -4.75e-07 3.555e-07 2.119e-07 -0.2819 0.7781 1
Subset.MCC.Mu.summary %>% select(-model_name) %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Extinction Estimates from Male-bias SS", split.table = Inf)
MCC BAMM Extinction Estimates from Male-bias SS
Parameter Estimate LCI UCI SE tval pval lambda
Sexual_selection_ppca 0.003441 -0.004092 0.01097 0.003844 0.8953 0.3707 1
log(range.size.m2) 0.0004368 -0.002251 0.003125 0.001371 0.3185 0.7501 1
bioclim4 -1.359e-05 -4.159e-05 1.441e-05 1.429e-05 -0.9511 0.3416 1
residuals.PC1 0.0001431 -0.002833 0.003119 0.001518 0.09423 0.9249 1
PC1.LIG -0.0009247 -0.003255 0.001406 0.001189 -0.7778 0.4368 1
NPP -2.11e-07 -1.169e-06 7.472e-07 4.889e-07 -0.4317 0.666 1

Table S15: HPD intervals are calculated for the above figure (i.e. models using Male-bias sexual selection measures). using this measure the models using \(\lambda_{DR}\) have model estimates that fall on the positive side 95 % of the time in 100 trees. For \(\lambda_{ND}\), there is a positive skew, however the 95 % HPD interval overlaps zero. \(\lambda_{BAMM}\) also shows a lesser positive skew. These intervals do not take into account the variance associated with each estimate and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter are in the same direction (+ or -). These intervals are calculated in the same way as in Table S7 and Table S9.

Subset.hpd.DR.top <- list()
Subset.DR.pgls.summary <- na.omit(Subset.DR.pgls.summary)
for(x in unique(Subset.DR.pgls.summary$Parameter)) {
Subset.hpd.DR.top[[x]] = hdi(Subset.DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.DR.top <- bind_rows(Subset.hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper"))


saveRDS(Subset.hpd.DR.top, 'data/Subset.hpd.DR.top.rds')

#For ND
Subset.hpd.ND.top <- list()
Subset.ND.pgls.summary <- na.omit(Subset.ND.pgls.summary)
for(x in unique(Subset.ND.pgls.summary$Parameter)) {
Subset.hpd.ND.top[[x]] = hdi(Subset.ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.ND.top <- bind_rows(Subset.hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Subset.hpd.ND.top, 'data/Subset.hpd.ND.top.rds')

###############
Subset.hpd.Lambda.top <- list()
Subset.Lambda.pgls.summary <- na.omit(Subset.Lambda.pgls.summary)
for(x in unique(Subset.Lambda.pgls.summary$Parameter)) {
Subset.hpd.Lambda.top[[x]] = hdi(Subset.Lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.Lambda.top <- bind_rows(Subset.hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) 
saveRDS(Subset.hpd.Lambda.top, 'data/Subset.hpd.Lambda.top.rds')

################
Subset.hpd.Mu.top <- list()
Subset.Mu.pgls.summary <- na.omit(Subset.Mu.pgls.summary)
for(x in unique(Subset.Mu.pgls.summary$Parameter)) {
Subset.hpd.Mu.top[[x]] = hdi(Subset.Mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.Mu.top <- bind_rows(Subset.hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Subset.hpd.Mu.top, 'data/Subset.hpd.Mu.top.rds')

Subset.hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "Subset DR HPD Interval")
Subset DR HPD Interval
  Sexual_selection_ppca log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower 0.00571 -0.0162 -0.000103 -0.000428 -0.00423 -2.58e-06
Upper 0.0578 -0.000597 2.21e-05 0.0125 0.00928 1.39e-06
Subset.hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "Subset ND HPD Interval")
Subset ND HPD Interval
  Sexual_selection_ppca log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.000504 -0.00042 -3.67e-06 -3.56e-05 -0.000212 -8.57e-08
Upper 0.00158 0.000122 8.27e-07 0.000453 0.000265 5.48e-08
Subset.hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "Subset BAMM Speciation HPD Interval")
Subset BAMM Speciation HPD Interval
  Sexual_selection_ppca log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.013 -0.0126 -9.65e-05 -0.00743 -0.00988 -2.09e-06
Upper 0.0309 0.0166 5.45e-05 0.019 0.00472 1.99e-06
Subset.hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "Subset BAMM Extinction HPD Interval")
Subset BAMM Extinction HPD Interval
  Sexual_selection_ppca log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.0516 -0.0222 -0.000424 -0.0341 -0.0185 -3.14e-06
Upper 0.0576 0.0306 0.00017 0.0195 0.0158 2.77e-06

Phylogenetic Path Analysis

We undertook a phylogenetic path analysis using a model of causal pathways a priori. To do this we used the phylopath function. Our path analysis was undertaken on the dataset subset to species with measures of male-biased sexual selection. Additionally, the path analysis was restricted to using one measure of speciation (\(\lambda_{DR}\)) and the MCC tree. Our reasoning for the causal links are described briefly below:

  1. \(\lambda_{DR}\) dependent upon Sexual Dichromatism: Sexual dichromatism would directly impact speciation if there is an association with the evolvability of plumage colour and speciation. This may more rapidly lead to niche divergence or reproductive isolation..
  2. \(\lambda_{DR}\) dependent upon Sexual Selection: This was a major hypothesis in our study (see introduction).
  3. \(\lambda_{DR}\) dependent upon Temperature Seasonality: More variable environments are broadly expected to impact the resources, niches and fitness of species; thus there is an expectation that more variable environments may impact speciation rate.
  4. \(\lambda_{DR}\) dependent upon Range Size: As previously seen in the PGLS models, range size shows a negative correlation with speciation rate. Expected reasons for this causal link include lower dispersal ability for species with smaller ranges correlating with lower gene flow and greater speciation. It may also reflect niche specialisation and thus, in line with the previous point, increase speciation and reduce gene flow.
  5. Sexual Dichromatism dependent upon Sexual Selection: Higher levels of sexual selection will promote showy males and drab females (sexual dichromatism).
  6. Sexual Dichromatism dependent upon Temperature Seasonality: Ecological aspects may facilitate sexual dimorphism as niche partitioning between the sexes causes sex-limited traits and colouration that are independent of mate choice or male-male competition.
  7. Sexual Selection dependent upon Temperature Seasonality: Greater environmental variability/stress may impact the usefulness (fitness effects) of sexual selection (see introduction).
  8. Range Size dependent upon Sexual Selection: Sexual selection may often depend on local adaptations and the alignment of sexually selected traits with traits increase local fitness. Thus, species with smaller ranges may have higher rates of local adaptation, with sexual selection further promoting the maintenance of local adaptations, restricted migration and thus smaller range sizes. Alternatively, if sexual selection has added fitness and evolvability benefits, it may enable species to colonise a greater range of environments leading to larger range size.
  9. Range Size dependent upon Temperature Seasonality: Anticipated to be more of a correlation than causal link; however a species range size may depend on the extent of environmental variation and availability of local resources.
#Set rownames to match tree
rownames(SS.subset) <- SS.subset$TipLabel
pruned.MCC.Subset.tree <- readRDS('data/pruned.MCC.Subset.tree.rds')
SS.subset$log.range.size <- log(SS.subset$range.size.m2)

SS.subset2 <- SS.subset %>% rename(
  DR = MCC.DR,
  SD = SDi,
  TS = bioclim4,
  RS = log.range.size,
  SS = Sexual_selection_ppca
)

models.Subset <- define_model_set(
  one = c(DR ~ SD, 
          DR ~ SS,
          DR ~ TS,
          DR ~ RS,
          SD ~ SS,
          SD ~ TS,
          SS ~ TS,
          RS ~ SS,
          RS ~ TS)
)

result <- phylo_path(models.Subset, data = SS.subset2, tree = pruned.MCC.Subset.tree, model = 'lambda')
  
#best_model <- best(result, boot = 500)
set.seed(1)
path.plot <- plot(x = best_model,
     type = "color",
     algorithm = 'gem',
     manual_layout = NULL,
     curvature = 0.1,
     colors = c("#b2182b", "#2166ac"),
     show.legend = F)

saveRDS(best_model, 'data/path_model.rds')
#Inspect the coefficients with their SE

# coef_plot(best_model, error_bar = "ci", reverse_order = TRUE) + 
#   ggplot2::coord_flip()+
#   ggplot2::theme_minimal()

# pdf("Figures/Path_Plot.pdf", width=8, height=8)
# path.plot
# dev.off()

# path.plot
best_model <- readRDS('data/path_model.rds')
coef_plot(best_model, error_bar = "ci", reverse_order = TRUE) + 
  ggplot2::coord_flip()+
  ggplot2::theme_minimal()

Figure S14: Path analysis standardised regression coefficients vary across relationships. Error bars are derived from confidence intervals through 500 bootstrapped iterations. Paths with error bars not overlapping zero are presented with an asterisks in Figure 3 within the manuscript.


Additional Figures and Tables

BAMM.df <- readRDS('data/BAMM.df.rds')
BAMM.Rates.100Trees <- do.call(rbind, BAMM.df)
BAMM.Rates.100Trees$mean.lambda <- BAMM.Rates.100Trees$mean.lambda
BAMM.Lambda.100Trees <- ggiraphExtra::summarySE(BAMM.Rates.100Trees, measurevar = "mean.lambda", groupvars = "TipLabel")

DR.list <- plyr::llply(es.list, function(x) {
  x %>% as.data.frame %>% rownames_to_column()
  })
DR.list <- do.call(rbind, DR.list)
DR.Lambda.100Trees <- ggiraphExtra::summarySE(DR.list, measurevar = ".", groupvars = "rowname")
DR.Lambda.100Trees <- DR.Lambda.100Trees %>% rename(DR.Lambda.100Trees, DR = ., TipLabel = rowname)

all.trees.join <- left_join(DR.Lambda.100Trees %>% select(TipLabel, DR), BAMM.Lambda.100Trees %>% select(mean.lambda, TipLabel), by = "TipLabel")

BAMM.Lambda.100Trees %>% ggplot(aes(x = DR.Lambda.100Trees$DR, y = mean.lambda, fill = mean.lambda))+
  geom_point(shape = 21)+
  geom_errorbar(aes(ymin = mean.lambda - 1.96*sd, ymax = mean.lambda + 1.96*sd), size = 0.02)+
  geom_errorbarh(aes(xmin = DR.Lambda.100Trees$DR - 1.96*DR.Lambda.100Trees$sd, xmax = DR.Lambda.100Trees$DR + 1.96*DR.Lambda.100Trees$sd), size = 0.02)+
  theme_minimal()+
  ylab("100 Trees Speciation Rate (BAMM)")+
  xlab("100 Trees Speciation Rate (DR)")+
  scale_fill_viridis_c()

# mean(DR.Lambda.100Trees$sd)
# mean(BAMM.Lambda.100Trees$sd)

Figure S15: Speciation Rate means from 100 trees using either the DR statistic (x axis) or BAMM (y axis). While there is a clear correlation (r = 0.75) there is variability with BAMM showing less heterogeneity. 95 % CIs are plotted for both axis from the 100 trees. Each point represents a species (n = 5,965).

restricted.data %>% ggplot(aes(x = MCC.DR, y = MCC.BAMM.model.data$mean.lambda, fill = MCC.DR))+
  geom_point(shape = 21)+
  theme_minimal()+
  geom_smooth(method = "lm")+
  ylab("MCC Speciation Rate (BAMM)")+
  xlab("MCC Speciation Rate (DR)")+
  scale_fill_viridis_c()

Figure S16: Similar to Figure S14, there is a correlation (r = 0.68) between \(\lambda_{DR}\) and \(\lambda_{BAMM}\) with BAMM results showing less heterogeneity. Each point represents a species (n = 5,965).

R Session information

sessionInfo() %>% pander()

R version 3.3.1 (2016-06-21)

**Platform:** x86_64-apple-darwin13.4.0 (64-bit)

locale: en_AU.UTF-8||en_AU.UTF-8||en_AU.UTF-8||C||en_AU.UTF-8||en_AU.UTF-8

attached base packages: parallel, grid, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: bindrcpp(v.0.2.2), ggtree(v.1.6.11), EBImage(v.4.16.0), phylopath(v.1.0.1.9000), stringr(v.1.2.0), ggraph(v.1.0.1.9999), stargazer(v.5.2.2), tibble(v.1.3.4), RColorBrewer(v.1.1-2), kableExtra(v.1.0.1), phangorn(v.2.4.0), HDInterval(v.0.2.0), MuMIn(v.1.42.1), coda(v.0.19-1), car(v.3.0-2), carData(v.3.0-1), BAMMtools(v.2.1.6), dplyr(v.0.7.6), ggExtra(v.0.8), reshape2(v.1.4.3), purrr(v.0.2.4), caper(v.1.0.1), mvtnorm(v.1.0-6), MASS(v.7.3-50), ggridges(v.0.5.0), brms(v.2.6.1), Rcpp(v.0.12.18), phytools(v.0.6-44), maps(v.3.2.0), gridExtra(v.2.3), geiger(v.2.0.6), lme4(v.1.1-18-1), Matrix(v.1.2-14), repmis(v.0.5), diversitree(v.0.9-10), ape(v.5.1), mgcv(v.1.8-24), nlme(v.3.1-128), tidyr(v.0.7.2), knitr(v.1.22.1), pander(v.0.6.2) and ggplot2(v.3.1.0.9000)

loaded via a namespace (and not attached): R.utils(v.2.7.0), tidyselect(v.0.2.4), htmlwidgets(v.1.2), combinat(v.0.0-8), munsell(v.0.5.0), animation(v.2.5), codetools(v.0.2-14), units(v.0.5-1), DT(v.0.4), ggiraphExtra(v.0.2.9.1), miniUI(v.0.1.1.1), withr(v.2.1.2), Brobdingnag(v.1.2-6), colorspace(v.1.3-2), uuid(v.0.1-2), highr(v.0.6), rstudioapi(v.0.7), stats4(v.3.3.1), officer(v.0.2.1), bayesplot(v.1.6.0), labeling(v.0.3), rstan(v.2.17.3), mnormt(v.1.5-5), bridgesampling(v.0.5-2), rprojroot(v.1.3-2), clusterGeneration(v.1.3.4), xfun(v.0.3), R6(v.2.2.2), markdown(v.0.8), locfit(v.1.5-9.1), ggiraph(v.0.4.2), bitops(v.1.0-6), assertthat(v.0.2.0), promises(v.1.0.1), scales(v.1.0.0), gtable(v.0.2.0), tidygraph(v.1.1.0), rlang(v.0.2.2), scatterplot3d(v.0.3-41), splines(v.3.3.1), lazyeval(v.0.2.1), mycor(v.0.1), broom(v.0.5.0), inline(v.0.3.15), prediction(v.0.2.0), yaml(v.2.2.0), abind(v.1.4-5), threejs(v.0.3.1), crosstalk(v.1.0.0), backports(v.1.1.2), httpuv(v.1.4.5), rsconnect(v.0.8.8), tools(v.3.3.1), psych(v.1.8.4), gplots(v.3.0.1), BiocGenerics(v.0.20.0), plyr(v.1.8.4), base64enc(v.0.1-3), viridis(v.0.5.1), deSolve(v.1.20), zoo(v.1.8-3), haven(v.1.1.2), ggrepel(v.0.8.0), magrittr(v.1.5), data.table(v.1.11.4), openxlsx(v.4.1.0), colourpicker(v.1.0), sjmisc(v.2.6.3), R.cache(v.0.13.0), matrixStats(v.0.54.0), hms(v.0.4.2), shinyjs(v.1.0), mime(v.0.5), evaluate(v.0.10.1), fftwtools(v.0.9-8), xtable(v.1.8-2), shinystan(v.2.4.0), rio(v.0.5.10), jpeg(v.0.1-8), readxl(v.1.0.0), rstantools(v.1.5.1), crayon(v.1.3.4), KernSmooth(v.2.23-15), minqa(v.1.2.4), R.oo(v.1.22.0), StanHeaders(v.2.17.2), htmltools(v.0.3.6), later(v.0.7.3), tiff(v.0.1-5), udunits2(v.0.13), expm(v.0.999-2), sjlabelled(v.1.0.5), tweenr(v.0.1.5), ppcor(v.1.1), subplex(v.1.4-1), readr(v.1.1.1), cli(v.1.0.0), quadprog(v.1.5-5), R.methodsS3(v.1.7.1), gdata(v.2.18.0), bindr(v.0.1.1), igraph(v.1.1.2), forcats(v.0.3.0), pkgconfig(v.2.0.2), numDeriv(v.2016.8-1), foreign(v.0.8-66), xml2(v.1.2.0), dygraphs(v.1.1.1.6), stringdist(v.0.9.5.1), webshot(v.0.5.0), rvg(v.0.1.7), rvest(v.0.3.2), snakecase(v.0.9.2), digest(v.0.6.16), msm(v.1.6.6), rmarkdown(v.1.10), cellranger(v.1.1.0), fastmatch(v.1.1-0), gdtools(v.0.1.6), curl(v.3.2), shiny(v.1.1.0), gtools(v.3.8.1), nloptr(v.1.0.4), jsonlite(v.1.5), viridisLite(v.0.3.0), lattice(v.0.20-33), loo(v.2.0.0), httr(v.1.3.1), plotrix(v.3.7-3), survival(v.2.39-4), glue(v.1.3.0), xts(v.0.11-0), zip(v.1.0.0), png(v.0.1-7), shinythemes(v.1.1.1), ggforce(v.0.1.1), stringi(v.1.1.6) and caTools(v.1.17.1)

References

Armenta, J. K., P. O. Dunn, and L. A. Whittingham. 2008. Quantifying avian sexual dichromatism: A comparison of methods. Journal of Experimental Biology 211:2423. Journal Article.

BirdLife International and Handbook of the Birds of the World. 2017. Bird species distribution maps of the world. http://datazone.birdlife.org/species/requestdis.

Dale, J., C. J. Dey, K. Delhey, B. Kempenaers, and M. Valcu. 2015. The effects of life history and sexual selection on male and female plumage colouration. Nature 527:367–370. Journal Article.

Del Hoyo, J., A. Elliott, and D. Christie. 2011. Handbook of the birds of the world (Vols. 8-16). Book, Lynx Edicions 2003-2011.

Harvey, M. G., G. F. Seeholzer, B. T. Smith, D. L. Rabosky, A. M. Cuervo, and R. T. Brumfield. 2017. Positive association between population genetic differentiation and speciation rates in new world birds. Proceedings of the National Academy of Sciences 114:6328–6333.

Jetz, W., G. H. Thomas, J. B. Joy, K. Hartmann, and A. O. Mooers. 2012. The global diversity of birds in space and time. Nature 491:444–448. Journal Article.

Quintero, I., and W. Jetz. 2018. Global elevational diversity and diversification of birds. Nature 555:246.

Rabosky, D. L., J. Chang, P. O. Title, P. F. Cowman, L. Sallan, M. Friedman, K. Kaschner, et al. 2018. An inverse latitudinal gradient in speciation rate for marine fishes. Nature 559:392.

Schliep, K. P. 2010. Phangorn: Phylogenetic analysis in r. Bioinformatics 27:592–593.

Title, P. O., and D. L. Rabosky. 2018. Diversification rates and phylogenies: What are we estimating, and how good are the estimates? bioRxiv 369124.